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ABSTRACT 

This  paper  is  a  report  on  a  continuing  effort  to  develop  a  physical-numerical 
model  suitable  for  use  in  predicting  large-scale  fields  of  low  cloudiness. 

The  need  to  incorporate,  within  such  a  model,  those  physical  processes  which 
are  characteristic  of  the  atmospheric  boundary  layer  has  been  noted  in  several 
recent  studies.  That  the  inclusion  of  such  processes  be  effected  without  imposing 
a  requirement  for  non-standard,  observational  data  has  been  adopted  as  a  restriction 
in  the  design  of  the  prediction  model.  A  further  design  specification  is  that  the 
model  will  be  used  in  conjunction  with  a  fine-mesh,  free-air  prediction  model. 

This  specification  relieves  us  from  the  need  to  consider  the  prediction  of  dynamical 
developments  within  the  bulk  of  the  atmosphere. 

Included  in  this  paper  are  a  derivation  of  a  basic  prediction  model  and  the 
results  of  a  test  series  of  12-hr  forecasts  made  with  various  versions  of  the  model 
using  synoptic  data  for  the  period;  12Z,  6  Feb. — OOZ,  7  Feb.  1964.  The  tests  were 
designed  to  assess  the  characteristics  of  the  model  and  to  indicate  areas  in  which 
more  study  is  required. 

The  model  tested  differs  broadly  from  previous  boundary  layer  models  in  its 
use  of  all  three  space  dimensions,  a  horizontal  space  mesh  of  150  km,  and  a  time 
step  of  15  min.  The  model  incorporates  the  computation  of:  eddy  fluxes  of  heat 
and  vapor,  the  transport  of  heat,  vapor,  etc.  by  ageostrophic  horizontal  winds,  the 
influence  of  terrain-  and  friction-induced  vertical  motion,  and  the  heat  and  mass 
exchanges  involved  in  water-substance  phase  changes. 

REVIEW  AND  APPROVAL 

Publication  of  this  technical  documentary  report  does  not  constitute  Air  Force 
approval  of  the  report’s  findings  or  conclusions.  It  is  published  only  for  the  exchange 
and  stimulation  of  ideas. 
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I 

INTRODUCTION 

This  report  documents  the  status  of  our  effort  to  develop  a  physical 
model  which  is  to  be  used  in  the  prediction  of  large-scale  fields  of  low  cloudi¬ 
ness.  A  previous  report  [1]  points  out  the  apparent  significance  of  those  physical 
processes  associated  with  the  structure  of  the  atmospheric  boundary  layer  for 
low-cloud  prediction.  Our  initial  assessment  of  the  feasibility  of  incorporating 
these  processes  in  a  physical-prediction  model  has  been  discussed  previously 
[17].  This  paper  reports  on  the  base-prediction  model  which  has  been  developed 
and  presents  the  results  of  a  series  of  tests  of  this  model  made  with  data  from 
one  synoptic  situation. 

It  should  be  differentiated  from  the  sub-synoptic  scale  boundary- layer 
models  discussed  in  subsection  2,  The  model  presented  here  has  been  designed 
to  be  used  in,  conjunction  with  a  free-air  prediction  model  developed  by  Air 
Weather  Service/  and  for  use  with  routine  meteorological  observations.  Our 
model  is  still  in  the  formative  stage  and  further  modification  and  testing  are 
planned  during  the  coming  year.  The  experimental  results  reported  here  are 
incompletely  documented  because  they  are  intended  to  indicate  only  the  general 
characteristics  of  the  model  and  areas  in  which  further  work  is  needed. 


We  refer  to  the  fine-mesh  model  under  development  by  Maj .  J.  Howcroft, 
1210th  Weather  Squadron,  Suitland,  Md. 
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II 


PREVIOUS  WORK 
1.  Cloud  Prediction  Techniques 

Other  tasks  involved  in  cloud  prediction  or  specification  under  Project  1 
of  the  433 L  Meteorological  Teclmique  Program  have  employed  an  empirical 
approach  to  the  problem  of  low-cloud  prediction.  Numerical-analysis  tech¬ 
niques,  based  on  synoptic  reasoning,  have  been  used  to  relate  the  occurrence, 
amount,  and  height  of  cloudiness  to  a  large  number  of  meteorological  quantities. 
A  relatively  recent  summary  of  this  work  has  been  reported  by  Cooley,  et  al.  . 
[6]. 

Another  approach  to  the  problem  is  based  on  the  numerical  integration  of  . 
simplified  conservation  equations  for  atmospheric  humidity  variables.  This 
approach  makes  use  of  air  trajectories  computed  from  the  wind  components 
which  are  forecasted  by  large-scale,  free-atmosphere,  dynamical  models 
[I,  22,  24j. 

Yet  another  approach  is  that  spelled  out  by  Smagor insky  [37]  in  which  an 
internally  consistent,  free  air,  dynamical  prediction  model  which  includes-  the 
condensation  process  is  employed  to  predict  the  entire  behavior  of  a  model 
atmosphere.  This  approach  was  initiated  by  Smagorinsky  and  Collins  [38]  some 
ten  years  ago.  The  methods  discussed  above  have  been  developed  largely 
because  it  did  not  prove  profitable  to  employ  the  full  dynamical  approach  in 
routine  forecasting  operations. 

All  three  techniques  discussed  above  have  a  common  weakness  —  their 
use  of  data  at  only  the  standard  pressure  levels.  It  can  be  seen  from  even  a 
cursory  examination  of  atmospheric  soundings  that  such  data  are  insufficient 
to  describe  the  distribution  of  humidity  in  the  atmosphere.  Further,  the  use 
in  these  techniques  of  geostrophic  or  balanced  wind  fields  to  compute  the  low 
level  transport  of  atmospheric  humidity  cannot  be  regarded  as  sound.  These 
factors,  together  with  the  neglect  of  other  boundary- layer  processes,  may  in 
part  explain  the  relative  lack  of  success  of  these  techniques  in  the  prediction 
of  low  cloudiness. 
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2.  Boundary-layer  Models 


Fisher  [14]  and  Estoque  [12],  in  continuation  of  some  initial  work  by- 
Pearce  [31],  approached  the  problem  of  numerical  prediction  of  the  local,  sea- 
breeze  circulation.  Their  treatment  of  this  phenomenon  led  naturally  to  rudi¬ 
mentary  modeling  of  certain  boundary  layer  processes.  Subsequently,  they 
embarked  on  studies  [13,  15]  of  boundary- layer  processes  without  the  complica¬ 
tion  of  local  circulations  in  the  wind  field. 

More  recently,  Pandolfo  [30,  31]  has  undertaken  the  formulation  of  a 
boundary-layer  prediction  model  for  use  in  predicting  short-period  changes 
in.  cloudiness,  visibility,  and  wind  velocity.  His  model  is  intended  for  use  with 
data  from  specially  designed  observation  networks. 

In  the  work  of  Estoque,  Fisher,  and  Pandolfo,  it  has  been  assumed  that 
the  boundary-layer  processes  are  confined  to  the  region  within  some  two 
kilometers  of  the  ground.  The  dependence  of  the  boundary-layer  phenomena 
upon  developments  in  the  free  air  above  that  level  is  accounted  for  by  ascribing 
the  task  of  providing  upper-boundary  conditions  for  the  boundary-layer  model 
to  an  independent,  free-air  prediction  model. 

Additionally,  Estoque  and  Pandolfo  employ  the  hypothesis  that  through  a 
shallow  layer  in  proximity  with  the  ground,  the  eddy  fluxes  of  momentum,  heat, 
and  vapor  are  constant  with  respect  to  height.  This  constant-flux  hypothesis 
is  discussed  in  the  book  by  Lumley  and  Panofsky  [25]. 

3.  Approach  Used  in.  Modeling 

The  low-cloud  forecasting  problem  being  studied  in  this  work  involves 
regarding  clouds  as  susceptible  to  adequate  specification  if  the  humidity  and 
temperature  fields  are  provided  on  a,  horizontal  grid-point  network  with  char¬ 
acteristic  spacing  of  .150  km.  This  horizontal-length  scale  and  the  characteristic, 
forecast-time  interval  of  one- half  day  differentiate  our  formulation  of  the  problem 
from  that  undertaken  in  the  work  discussed  in  Subsection  2.  We  do  model  the 
vertical  structure  of  the  lower  atmosphere  in  the  fashion  of  Estoque  and  Pandolfo. 


This  is  done  to  permit  adequate  vertical  resolution  for  the  computation  of  the 
boundary-layer  process’' of  eddy  diffusion.  Such  emphasis  on  boundary-layer 
processes  has  not  been  given  in  the  other  techniques  discussed  in  subsection  1. 

Our  model  has  been  designed  to  require  only  routine  observational  data  for 
its  implementation.  Because  of  this  design  requirement,  we  did  not  choose  to 
attempt  the  integration  of  the  time-dependent,  horizontal-momentum  equation. 
The  horizontal-wind  components  are  computed  from  the  diagnostic  equations, 
which  result  if  one  assumes  an  instantaneous  adjustment  of  the  wind  to'  create 
a  balance  among  the  pressure  gradient,  Coriolis,  and  frictional  forces. 

Consequent  upon  this  scheme  for  computing  the  wind  field  is  our  use  of 

empirical  relationships  between  the  geostrophic  wind  and  the  friction  velocity. 

Additionally,  it  was  necessary  to  use  simple  formulations  for  the  distribution  of 

2 

the  mixing  coefficient  throughout  the  transition  layer. 

The  large-scale  character  of  the  model  has  also  required  our  omission 
of  the  interface  energy  balance  technique  [13]  for  computing  the  surface 
temperature.  We  propose  instead  to  employ  an  empirical  approach  to  the 
specification  of  this  boundary  value. 


The  boundary  layer  is  divided  into  a  layer  of  constant  flux  (contact  layer)  and 
a  transition  layer. 
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III 

TRANSITION- LAYER  PREDICTION  EQUATIONS 


4.  Introduction 

The  coordinate  system  used  for  expressing  the  model’s  governing  equations 
is  a  modification  of  the  quasi-Cartesian  system  [19].  Because  the  pertinent 
physical  boundary  is  the  surface  of  the  ground  or  ocean,  and  not  a  fictitious  mean 
sea  level  we  choose  as  the  vertical  coordinate  a  parameter,  z,  which  takes  on  the 
value  zero  at  the  ground  level.  Thus,  if  with  value  zero  at  mean  sea  level, 
is  the  customary  vertical  coordinate  of  the  quasi-Cartesian  system  and  if  E(x,  y) 
is  a  function,  giving  the  elevation  of  the  ground  above  mean  sea  level  at  a  point 
with  horizontal  coordinates,  (x, y),  then  z  is  given  in  terms  of  |  and  E  by 

z  =  ^  -  E(x,y).  (III-l) 

A  detailed  derivation  of  the  form  of  the  various  differential  operators  in  the 
modified  coordinate  system  is  given  in  Appendix  I, (note  the  difference  in  the 
definition  of  E) , 

In  accordance  with  the  assumptions  regarding  the  structure  of  the  boundary 
layer  outlined  earlier,  the  top  and  bottom  of  the  transition  layer  are  the  coordinate 
surfaces, 

z  =  H  (2  km)  and  z  =  h(50m),  (111-2) 

respectively.  In  the  formulation  of  the  prediction  equations,  we  will  call  for  the 
specification  of  boundary  conditions  on  these  coordinate  surfaces.  In  the  case 
of  the  upper  surface,  the  boundary  conditions  will  be  taken  as  prescribed  by 
means  of  predictions  made  with  an  independent,  free  air  model.  This  is  a 
temporary  solution  to  the  ultimate  problem  of  designing  an  internally  consistent, 
physical  model  for  the  entire  troposphere.  The  required  lower  boundary  condi¬ 
tions  will  be  obtained  through  the  use  of  the  equations  applicable  to  the  contact 
layer. 

Additional  boundary  conditions  will  be  required  on  the  open  lateral  boundaries 
of  the  region  to  which  the  model  is  applied.  These  conditions  require  specifica¬ 
tion  of  the  horizontal  advection  of  the  various  quantities  into  the  region.  This 
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information  cannot  be  given  in  realistic  applications  of  the  model  to  s3moptic 
data;  thus,  one  is  forced ^to  provide  some  approximation  for  these  quantities 
and  to  accept  the  inward  propagation  of  error.  The  propagation  speed  of  such 
errors  will  be  the  speed  of  the  local  wind,  and  the  errors  will  not  rapidly  cor¬ 
rupt  the  forecasts  in  the  interior  of  a  sufficiently  large  region. 

The  wind  velocity  and  mixing  coefficients  appear  in  the  equations  derived 
below.  These  quantities  will  be  determined  through  a  series  of  relations  which 
are  both  theoretical  and  empirical.  They  will  not  be  governed  by  simple  pre¬ 
diction  equations. 


5.  Wind  Specification 

The  horizontal  wind  coniponents,  u  and  v,  are  taken  to  be  the  solutions  of 
the  following  simplified  equations  of  motion: 


_f 

K 


(V 


(III-3a) 


K 


(u 


u  ) , 
g 


(III-3b) 


where  f  is  the  coriolis  parameter  and  K  is  the  momentum  mixing  coefficient 
which  is  assumed  to  be  independent  of  z.  The  geostrophic  wind  components, 
u  and  V  ,  are  assumed  to  be  linear  functions  of  height. 


H  (H  -  z)  ,  H  0 
u  =  u  -  -  (U  -  U  )  , 

g  g  H  '  g  g'’ 


(III-4a) 


H  (H  -  z)  H  0^ 
g  g  H  '  g  g' 


(III-4b) 


H  H  0  0 

where  u  and  v  are  the  geostrophic  wind  components  at  z  =  H,  and  u  and  v 
g  g  '  g  g 

are  the  geostrophic  wind  components  at  z  =  0. 

The  solution  of  Eqs,  III-3a  and  3b  is  specified  by  requiring  that 


u  u 


g 


g 


as 


as 


(III-5a) 

(in- 5b) 
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and  that 


u  =  U  (X;  y,  t)  =  U 

at 

z.  =  h  , 

(Ill- 6a) 

V  =  V(x,y,tl  -  V 

at 

z  =  h  . 

(III-6b) 

One  may  verify  that  the  appropriate  solution  to  the  given  system  is 
u  =  u  -t-  e  {(U  -  )  cos  [Oi  {7  -  h)|  +  (V  -  v^)  sin.  (o;  (z  -  h)]} 

s  s  s 


(III-7a) 


V  =  V  +  e  {(V  -  v^i  cos  [O!  (z  -  h)]  -  (U  -  u^)  sin  [o^  (z  -  h)]) 

s  s  s 


(III-7b) 


where  u^  and  are  the  values  of  u  and  v  at  z  ~  h,  and  o'  is  given  by 
g  g  g  g 


a  =  — 
2K 


(III-8) 


The  boundary  conditions  at  z  =  h  are  prescribed  through  the  appropriate 

contact  layer  formulas  given  later.  If  one  were  to  take  u  and  v  as  constants, 

g  g 

and  U  -  V  =  0  at  z  =  h  =  0.  the  solution  would  reduce  to  an  Ekman  spiral.  The 
simple  linear  dependence  of  the  geostrophic  wind  upon  height  allows  the  super¬ 
position  of  a  thermal  shear  upon  the  spiral.  The  specification  of  boundary 
conditions,  U  and  V,  at  z  =•  h,  imposes  a  more  realistic  shear  between  the  ground 
and  z  =  h  than  would  result  from  an  Ekman  spiral  solution.  This  may  be  attributed 
to  the  marked  variation  of  the  mixing  coefficient  with  height  in  the  contact  layer. 
The  value  of  the  mixing  coefficient  used  in  evaluating  the  wind  components  varies 
with  horizontal  position  and  time,  in  accordance  with  the  formulas  subsequently 
derived  for  the  contact  layer. 

The  vertical  wind  in  the  transition  layer  has  been  divided  into  two  components 
in  accordance  with  the  structure  of  the  coordinate  system  being  employed  (see 
Appendix  I).  A.  frictionally-induced  component  is  obtained  from  mass  continuity 
considerations  if  we  neglect  the  individual  change  in  particle  density  (a  reasonable 
approximation  (21])  in  the  continuity  equation. 

This  wind  component  is  assumed  to  vanish  at  z  =  h  and  is,  therefore,  given 
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by  the  relation 


w 


^  ^  ^  9v : 

f  I  ax  dy  I 
h 


(III-9) 


where  u  and  v  are  the  horizontal  wind  components  derived  earlier.  Although  the 
horizontal  winds  are  distributed  in  a  more  complex  fashion  than  that  given  by  an 
Ekman  spiral  solution,  one  may  still  qualitatively  associate  upward  currents 
with  low  pressure  centers  and  downward  flow  with  high  pressure  centers. 

The  second  component  of  the  vertical-wind  velocity  is  related  to  the  varia¬ 
tion  of  the  ground  elevation  from  mean  sea  level.  We  define  the  terrain-induced 
velocity,  w,  as, 

^  aE  aE 

w  =  u—  +  V  —  ,  (III-IO) 

ax  ay  '  . 

where  u  and  v  are  functions  of  x,y,  z,  and  t,  which  are  given  by  Eqs.  III-7a  and  7b 
This  velocity  will  vary  with  height  due  to  the  height  dependence  of  u  and  v. 
Because  of  the  coordinate  system  being  employed,  this  component  of  the  wind 
will  appear  explicitly  only  in  the  thermodynamic  energy  equation,  in  which  it  is 
associated  with  adiabatic  temperature  changes. 


6.  The  Specification  of  Mixing  Coefficients 

By  definition,  the  transition  layer  is  a  region  in  which  the  direct  influence 
of  the  ground  on  atmospheric  properties  gradually  diminishes.  This  layer  can 
be  taken  to  be  of  fixed  depth  only  with  respect  to  certain  scales  of  atmospheric 
phenomena.  We  regard  the  large-scale,  vertical  eddy  exchange  as  being  pre¬ 
dominantly  confined  to  a  comparatively  thin  layer  above  the  ground. 

The  influence  of  vertical  eddy  exchange  may  be  expressed  in  a  prediction 
equation  through  a  term  of  the  form. 


9z  9z  ’ 


(III-ll) 


where  9  represents  an  arbitrary  atmospheric  property,  and  K  is  the  mixing 
coefficient  appropriate  to  that  property.  The  product  of  K  and  the  derivative  of 
0  represents  the  eddy  flux  of  the  property. 

This  term  may  be  expanded  by  carrying  out  the  indicated  differentiation; 


—  IK  —1  =  —  —  +  K  — 
9z  9z  9z  9z  2’ 

9z 


(III-12) 
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The  second  term  on  the  right  hand  side  of  this  equation  is  of  a  form  encountered 
in  the  molecular  theory  of  heat  transfer.  The  action  of  K  in  this  term  is  restricted 
to  reducing  (increasing)  local  maxima  (minima)  in  the  vertical  distribution  of  0  at 
a  rate  dependent  upon  the  magniitude  of  K.  In  the  first  term  on  the  right  hand 
side  of  Eq.  111-12,  we  see  a  form  suggestive  of  that  associated  with  convective 
transfer  of  the  property  0  by  non-eddy  vertical  air  currents.  Suppose,  for 
concreteness,  that  0  is  a  property  which  increases  with  z.  Then,  if  9K/9z  is 
negative,  this  term  will  contribute  to  a  local  decrease  in  0  in  the  same  way  as 
would  an  upward  wind  current. 

If  one  reflects  on  the  mixing-length  concept  of  turbulent  exchange  [20],  it 
seems  reasonable  to  anticipate  that  the  coefficient  K  would  tend  toward  a  limiting 
value  with  increasing  distance  from  the  ground.  The  size  of  this  value  is  expected 
to  depend  upon  the  two  factors,  intensity  of  the  turbulent  motion  induced  at  the 
ground,  and  the  distribution  of  the  property  being  transferred.  Additionally,  the 
possibility  of  the  atmospheric  static-stability  structure  being  such  as  to  inhibit 
the  vertical  extent  of  ground-produced  eddy  motion  can  be  expected  to  influence 
the  vertical  distribution  of  K. 

We  do  not  possess  an  adequate  theory  for  prescribing  the  precise  behavior 
of  the  mixing  coefficient  in  the  transition  layer.  Nonetheless,  we  have  attempted 
a  formulation  of  its  variation  based  on  the  ideas  outlined  above. 

To  begin  with,  a  basic  value  of  the  coefficient  is  prescribed  at  the  base  of 
the  transition  layer  in  accordance  with  the  value  derived  from  the  contact-layer 
formulas.  This  value  gives  an  estimate  of  the  intensity  of  the  mixing  process 
active  in  the  contact  layer,  and  has  been  used  to  specify  in  two  ways  the  coefficient, 
K,  throughout  the  remainder  of  the  transition  layer. 

First,  it  was  assumed  that  this  value  remains  unchanged.  It  is  this  assumption 
that  was  used  in  the  simplified  equations  of  horizontal  motion.  The  second  assump¬ 
tion  allowed  for  a  linear  decrease  of  this  coefficient  with  height  to  a  residual  value 

3 

of  10  cgs  units  at  the  top  of  the  transition  layer.  In  order  to  allow  for  the  influence 
of  atmospheric  stability,  we  introduced  a  weak  dependence  upon  thermal  stability 
from  500  m  above  the  ground  to  the  top  of  the  transition  layer.  The  restriction 
to  elevations  above  500  m  was  prompted  by  the  results  of  some  initial  experiments. 
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Using  idealized  data  with  a  two-connponent,  sinusoidal,  diurnal  ground-temper¬ 
ature  wave,  and  permitting  the  stability  dependence  to  extend  down  to  the  base 
of  the  transition  layer  resulted  in  considerable  distortion  of  the  diurnal  tempera¬ 
ture  wave.  This  second  assumption  on  the  mixing  coefficient  was  not  introduced 
into  the  horizontal  wind  computation  because  of  the  difficulty  and  computation 
time  involved  in  computing  a  numerical  solution  of  the  governing  equations. 

Denoting  the  value  of  the  mixing  coefficient  at  the  top  of  the  contact  layer 
by  K,  the  second  assumption  on  the  distribution  of  K  may  be  expressed  as, 

K  =  10^  +  K  1(H  -  z)/(H  -  h)]  ,  h  <  z  <  500  m,  (III-13a) 

3  +4  9T 

K  =  10  +  K  1(H  -  z)/(H  -  h)]  (0.88  -  0.20  -10  — )  ,  500m  <  z  <  H. 

(III-13b) 
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We  also  imposed  an  upper  limit  of  10  cgs  units  on  the  value  of  K, 


7.  The  heat  Transfer  Equation 

If  the  only  diabatic  energy  source  is  the  convergence  of  the  vertical-eddy 
heat  flux,  the  first  law  of  thermodynamics  and  the  equation  of  the  state  for  an 
ideal  gas  may  be  combined  to  yield  the  equation. 


30  30  30  30  3  30 

—  —  -  u —  -  V —  -  w —  +  —  K  — 
3t  3x  3y  3z  3z  H  3z 


(III-14) 


where  0  is  the  potential  temperature  and  K  is  the  mixing  coefficient  for  sensible 

H 

heat.  The  potential  temperature  is  related  to  the  temperature,  T,  and  the  pres¬ 
sure,  p,  through  the  expression; 

-K 

0  =  T  (^  ,  (III-15) 


where  P  is  a  standard  pressure  and  k  is  the  ratio -of  the  gas  constant  for  dry  air, 
R,  to  the  specific  heat  at  constant  pressure  for  dry  air,  c^.  The  hydrostatic 
equation. 


3z 


gE_ 

RT’ 


(III-16) 


where  g  is  the  acceleration  of  gravity,  may  be  used  to  rewrite  Eq.  III-14  as 


dt 


KT 

^  +  -9 

K 

^  ^  _g  ■ 

+  •  K 

^  _g_' 

p 

dt  3z 

H 

3z  c 

P  J 

c  T  H 

P 

3z  c 

L  p 

(III-17) 
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The  first  term  on  the  right  can  be  approximated  (see  Appendix  I)  by  the  relation, 


.<cT  dp 
p  dt 


XT'. 

- (W  +  Wl 

p  ■ 


dz 


(w  +  w) . 


(III-18) 


The  last  term  on  the  right  of  Eq.  111-17  is  negligibly  small  except  in  cases  of 
extreme  instability,  which  we  will  not  expect  to  encounter  in  application  of  the 
model.  By  introducing  these  approximations,  Eq.  111-17  may  be  put  into  the 
form , 


91  _  _  .9T 

9t  ^  9x 


9T 

V  -  W 

9y 


9T 

—  + 

9z  c„ 


w  +  K 


H  9z  c_ 


•  (III-19) 


In  addition  to  the  diabatic  influence  of  the  eddy  heat  flux,  we  have  incorporated 
into  Eq.  111-19  a  term,  (dT/dt)w,  representing  the  heat  exchange  associated 


with  water  substance  phase  changes.  The  procedure  for  computing  this  diabatic 
effect  is  discussed  in.  Subsection  9. 


8.  Water  Substance  Prediction 

We  will  consider  water  to  be  present  in  the  atmosphere  in  two  phases: 
liquid  and  vapor.  The  ratio  of  the  mass  density  of  water  vapor  to  the  total  mass 
density  of  moist  air  (vapor,  dry  air,  and  liquid  water)  is  defined  as  the  specific 
humidity  and  denoted  by  the  letter,  q.  We  will  simila.rly  define  a  quantity,  r,  as 
the  ratio  of  the  mass  density  of  water  substance  (vapor  plus  liquid)  to  the  total 
mass  density  of  moist  air.  This  quantity  will  be  referred  to  as  the  specific 
moisture. 

If  we  imagine  a  parcel  of  moist  air  to  move  in  space,  then  the  specific 
moisture  will  be  invariant  following  such  a  movement  under  the  following 
conditions; 


(a)  no  molecular  or  eddy  diffusion  occurs,  and 

(b)  no  precipitation  enters  or  leaves  the  parcel. 

The  water  may  change  phase  during  the  motion  of  the  parcel,  but  the  total 
moisture  will  be  invariant.  Thus,  the  equation  governing  r  may  be  written 

^=0.  (III-19) 
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The  specific  humidity  of  the  parcel  may,  of  course,  change  under  these 
conditions  in  response  to'i water-phase  variations.  If  we  denote  the  depletion 
rate  of  water  vapor  in  response  to  condensation  by,  (dq/dt)^  we  may  write  the 
following  equation  for  specific  humidity; 


dt 


dt 


(III- 20) 


Because  the  atmospheric  motion  is  turbulent,  eddy  diffusion  of  water  sub¬ 
stance  will  occur.  It  will  be  assumed  that  the  eddy  niixing  coefficients  for  both 
vapor  and  liquid  are  equal  and  that  the  liquid  and  vapor  follow  the  air  motion 
[27],  Equations  III-19  and  III- 20  then  become,  upon  expansion  of  the  particle 
derivatives  and  the  introduction  of  the  eddy  exchange  mechanisms. 


3q  9q  3q  8q 

=  -  U —  -  V~^  -  W~^  +  .  -- 

3t  3x  3y  3z  3z  1  w  3z 


^1- 


[dtj 


(III- 21) 


and 


3r  3r  3r  3r  3  3r 

—  =  -  u —  -  V —  -  w —  +  —  K  — 

3t  3x  3y  3z  3z  w  3z 


(III- 22) 


The  coefficients  in  Eqs.  III-21  and  22,  and  in  Eq.  III-19  are  com¬ 
puted  from  Eqs.  III-13a  and  13b.  This  procedure  is  followed  because  of  its 
simplicity  and  the  unavailability  of  any  theoretically  sound  alternative  within 
the  framework  of  this  model. 


9.  Computation  of  Water  Phase  Changes  and  Diabatic  Heating 

The  equation  governing  specific  humidity  contains  a  term,  (dq/dt)  ,  which 

c 

represents  the  source  or  sink  of  water  vapor  associated  with  phase  changes  in 
the  water  content  of  the  air  parcel.  Similarly,  the  heat  release  associated  with 
these  phase  changes  appears  in  the  thermodynamic  energy  equation  as  the  diabatic 
source  term  (dT/dt)^.  In  order  to  compute  these  two  terms,  we  must  introduce 
approximating  assumptions.  The  scheme  employed  is  based  on  the  method  used 
by  Fisher  [15]  and  modified  by  McDonald  [27]. 

If  the  parameters,  q,  r,  and  T,  and  the  other  model  variables  are  known 
at  a  given  time,  the  equations  governing  q,  r,  and  T,  omitting  the  source  terms 


(dq/dt)  and  (dT/dt)  ,  may  be  used  to  obtain  a  first  approximation  of  the  values 
0  0  ^  0 

q  ,  r  ,  and  T  at  the  next  time  step.  The  saturation  value  of  specific  humidity 

at  a  given  point  can  be  regarded  as  a  function  of  the  temperature  alone  if  the 

pressure  dependence  is  appropriately  approximated.  Thus,  with  the  value  T*^ 

we  may  associate  a  saturation  value  of  specific  humidity,  q^.  An  analysis  of 

0  0  0  0  i 

the  relative  values  of  the  four  quantities  —  q  ,  r  ,  T  and  q  —  must  be  made 

s 

to  decide  if  the  values  should  be  adjusted  to  account  for  the  occurrence  of  phase 
changes  and  associated  diabatic  temperature  changes.  Before  discussing  the  logic 
of  the  adjustment  procedure,  we  first  derive  the  equations  for  evaluating  the 
approximate  value  of  the  adjustment. 

For  a  change  of  water  between  the  liquid  and  vapor  phases,  the  Clausius- 
Clapeyron  equation  may  be  written  to  good  approximation  [33]  as 


d(lne^)  =  m^Ld(T)/(R*T  ), 


(III-23) 


where  e^  is  the  saturation  pressure  of  the  water  vapor  in  equilibrium  with  a 

plane  water  surface  at  a  temperature  of  T  degrees  absolute,  m^  is  the  molecular 

weight  of  water  vapor  (18),  L  is  the  latent  heat  of  vaporization  (cal  gm  ^),  and  R* 

is  the  universal  gas  constant  (1.986  cal  deg  ^).  The  latent  heat,  L,  is  a  function 

of  the  vapor  temperature.  This  equation  may  be  used  to  derive  an  approximate 

relationship  between  the  saturation  specific  humidity,  q  ,  and  the  temperature. 

s 

The  value  of  q  is  closely  given  by 
s 


m  e 
w  s 

q  = - 

s  m^  p 


(in- 24) 


where  m^  is  the  molecular  weight  (-  29)  of  the  gas  mixture  referred  to  as 
“pure  dry  air”,  and  p  is  the  partial  pressure  of  this  mixture  in  the  parcel.  It 
follows  from  Eqs.  III-23  and  24  that 


d  (In  e  )  =  d  (In  q  )  +  d  (In  p)  =  m  L  d(T)/ (R*T  )  , 
s  s  w 


(III- 25) 


which,  for  normal  atmospheric  conditions  may  be  written  to  good  approximation 


m  Lq 

,  w  s  _ 

dq  =  - dT. 

R*T^ 


(III- 26) 
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Equation  III- 26  is  one  of  the  two  equations  which  will  be  employed  in  the  satura¬ 
tion  adjustment  computation. 

For  an  adiabatic-isobaric,  water  phase  change  (liquid  - — ►  vapor)  in  a  moist 
air  parcel,  the  first  law  of  thermodynamics  may  be  written  in  approximate  form 
as 


c  dT  =  -Ldq 
P 


(III-27) 


where  c^  is  the  specific  heat  at  constant  pressure  of  dry  air  and  L  is  the  latent 
heat  of  vaporization.  This  is  the  second  of  the  two  equations  which  are  required 
for  the  computation  of  the  saturation  adjustment. 

We  may  now  indicate  the  computational  form  of  these  two  equations.  Recalling 
that  a  superscript,  0,  implies  that  the  symbol  stands  for  the  first  approximation 
to  the  value  of  the  parameter  (obtained  as  outlined  earlier),  and  using  an  unscripted 
symbol  to  stand  for  the  value  of  the  parameter  after  the  saturation  adjustment" 
has  been  applied  to  it,  we  may  write  Eq.  III-26  as 


0  ^ 
'  (T»)^ 


(T  -  T  ) 


(III-28) 


and  Eq.  III-27  as 


b(T  -  T®)  =  -  (q  -  q®), 

where  a  and  b  are  the  quantities  (m^^  and  (Cp/L)  respectively. 


(III-29) 


Equation  III- 28  may  be  regarded  as  a  straight-line,  tangent  approximation 
0  0 

at  the  point  (q  ,  T  )  to  the  curve  representing  the  solution  of  the  Clausius- 
s 

Clapeyron  equation  in  a  (q,  T)  plane.  The  parameters  a  and  b  will  be  taken  to 
be  constants;  consequently,  Eq.  III-29  is  a  linear  approximation  to  the  general 
curvilinear  solution  to  Eq.  III-27. 

One  may  now  regard  Eqs.  Ill- 27  and  28  as  a  pair  of  simulianeous  linear  -  - 
equations  in  the  two  unknowns,  T  and  q.  Using  Eq.  III-29  in  Eq.  III-28,  one  may 


obtain. 


T^  + 


0  0 
q 


(III-30) 
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Similarly,  using  Eq.  Ill- 28  in  Eq.  Ill- 29,  one  finds  that 


0  0 
q  -  q^ 


q  =  q  -  b 


(III-31) 


aq^ 


+  b 


The  parameters  a  and  b  have  the  values 

a  -  5.393  X  10^  (deg), 

4  -1 

b  =  4.017  X  10  (deg  , ). 


(III-32) 

(III-33) 


The  logic  for  appl3nng  the  adjustment  Eqs.  III-30  and  31  will  now  be  presented 
using  the  symbols 


and 


,  0  0,^  2 

,  _ 

T  0  ^^0  2 

aq  +  b  T  . 
s 


c  =  -  b 
q  T 


(III-34) 


(III-35) 


The  first  step  is  to  compare  q^  and  q^ .  If  q^  is  equal  to  q^  ,  no  adjustment 

S  .  •  s 

is  required  and  one  simply  writes 

T  =  T^ 

0 


q  =  q 

0 

r  =  r  . 


(III-36) 


If  q^  is  greater  than  q^  ,  condensation  should  take  place  and  an  examination  of  the 

availability  of  water  vapor  must  be  included.  If  q^  is  greater  than  |c  1,  one 

Q 

writes 


q  =  q  +  c^ 


T  =  T  +  c. 


0 

r  =  r  . 


(III-37) 
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Note  in  this  case  that  c  <0  and  >  0,  and  that  ic  1  refers  to  the  absolute 
value  of  c^.  If  however,  jc^i  is  greater  than  or  equal  to  q^,  one  writes  the 
adjusted  values  as 

q  =  0 

r  =  r^  (III-38) 

T  -  T®  +  (q°/b)  . 


The  alternative  which  leads  to  Eq.  III-38  is  not  likely  to  occur,  as  it 
implies  the  complete  condensation  of  all  the  vapor  in  the  air.  We  include  it 
here  for  logical  completeness  alone. 

Finally,  if  q^  is  less  than  q^,  one  may  expect  evaporation  of  the  liquid  water 
in  the  air  parcel.  If  (r^  -  q^)  is  zero  (or  negative  for  programming  purposes) 
there  is  no  liquid  water  in  the  air  parcel  and  consequently  no  saturation  adjust-  . 
ment  is  called  for.  One  simply  writes  Eq.  III-36  again.  But  if  (r^  -  q^)  is 
positive,  there  is  liquid  water  available  for  evaporation.  (Note  that  in  this  case 
c^  >  0  and  c,^,  <  0).  One  must  then  inquire  if  c^  is  greater  than  (r*^  -  q*^); 
if  it  is  not,  one  writes 


T 


+  c 


T 


c 

q 


(III-39) 


but  if  it  is  greater,  then  one  can  only  permit  the  evaporation  of  the  liquid  water 
which  is  available.  Consequently,  the  adjustment  equations  are 


T  = 


q^/b . 


(III-40) 


The  examination  of  these  alternatives  exhausts  the  possible  situations  and 
leads  finally  to  the  specification  of  values  of  temperature  and  specific  humidity 
that  have  been  appropriately  modified  for  the  possibility  of  water  phase  changes. 
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T-r-^TT-r^ — n  nWF"l'f''T 


As  a  final  point  in  this  section,  we  indicate  the  formula  used  to  evaluate 

from  T^.  In  order  to  avoid  the  necessity  for  integrating  the  hydrostatic 
s 

equation  for  the  sole  purpose  of  evaluating  this  relationship,  we  have  assumed 
that  the  pressure  at  a  point  is  an  invariant  function  of  the  elevation  of  that 
point  above  mean  sea  level.  Using  the  ICAO  standard  atmosphere,  we  found 
that  this  dependence  could  be  written  to  good  approximation  as. 


P(x.y.  z) 


(1013  -  1.065  X  10 


(z  +  E(x,y)}  [mb]. 


(III-41) 


in  where  z  is  the  elevation  in  cm  of  the  point  above  the  ground  and  E  is  the 
elevation  in  cm  of  the  ground  above  mean  sea  level.  Then,  using  Tetens  approxi¬ 
mating  formula  for  saturation  vapor  pressure  quoted  in  [20],  and  Eq.  III-24,  we 
obtained  the  formula. 


3.8  X  10 


-3 


-6 


[1.013  -  1.065  X  10  (z  +  E)] 


exp  / 17.25 


T®  - 

1 

273 

T^  - 

35.7 

(III-42) 


10.  Cloud  Specification 

The  prediction  equations  derived  above  provide  the  basis  for  forecasting 
the  distribution  of  the  relative  humidity  and  liquid  water  content  of  the  air. 
Although  one  might  wish  to  assign  to  these  predicted  quantities  a  unique  inter¬ 
pretation  in  terms  of  cloud  occurrence,  such  a  procedure  is  not  apt  to  give  very 
accurate  results. 

In  view  of  the  experience  obtained  in  large-scale  dynamical  and  statistical 

cloud  prediction  studies  [1,  6,  37]  it  seems  likely  that  superior  results  will  be 

obtained  by  developing  a  probabilistic  interpretation  of  the  forecast  parameters 

in  terms  of  cloud  occurrence.  In  Smagorinsky’s  paper,  the  problem  of  specifying 

a  cloud  amount— relative  humidity  relationship  is  discussed.  The  results  of  a 

limited  study  indicate  that  a  linear  relation  between  these  parameters  could 

explain  a  significant  proportion  of  the  observations.  For  low  clouds,  cloud  amount 

in  tenths,  denoted  by  c,  was  related  to  the  ratio  of  q  to  q  by  the  linear  expression, 

s 

c  =  3.25  (q/q  )  -  1.95,  (q/q  )  ^  .6 .  (III-43) 

s  s 
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The  problem  of  specifying  this  relationship  occurs  not  only  in  the  interpretation 
of  forecasts,  but  also  in  the  analysis  of  the  initial  data.  Continuing  research 
into  these  analyses  and  forecast-interpretation  problems  will  eventually  lead 
to  an  adequate  solution,  but  an  account  of  these  efforts  is  beyond  the  scope  of 
this  paper.  We  have  employed  techniques  which  are  in  general  agreement  with 
Eq.  III-43  in  discussing  forecasts  made  with  this  model  and  in  performing  the 
analysis  of  initial  data  for  these  forecasts. 

11.  Radiative-heat  Flux  and  Precipitation 

Two  significant  physical  processes  were  omitted  in  formulating  the 
equations  governing  the  atmosphere  in  the  transition  layer.  It  is  a  priori 
desirable  to  incorporate  a  treatment  of  radiative -heat  transfer  and  precipitation 
into  a  physical  model  of  the  atmospheric  boundary  layer.  The  neglect  of  these 
processes  within  our  model  is  therefore  deserving  of  some  comment. 

We  shall  restrict  our  discussion  in  this  section  to  the  radiative-heat 
transfer  in  the  transition  layer.  The  topic  will  be  discussed  separately  as  it 
affects  our  treatment  of  the  contact-layer  equations. 

The  boundary-layer  models  developed  by  Estoque  [11]  and  Pandolfo,  Cooley, 
and  Atwater  [32]  incorporate  procedures  for  computing  the  convergence  of  the 
infrared  radiative  flux.  Neither  model  has  as  yet  been  tested  on  a  case  in 
which  clouds  were  present  within  the  boundary  layer.  The  basic  treatment  of  the 
radiation  process  follows  the  work  of  Elsasser  [10,11].  In  Estoque’s  model,  a 
numerical  procedure  developed  by  Brooks  [4]  and  adapted  by  Elliott  and  Stevens 
[9]  is  used.  In  Pandolfo’s  model,  the  Elsasser  tables  [11]  are  stored  in  the 
computer,  and  the  flux  computed  by  a  numerical  integration  analogue  of  the  area¬ 
measuring  graphical  technique. 

It  must  be  emphasized  that  both  of  these  sub-synoptic  scale  models  are 
oriented  toward  single-station  forecasting.  This  permits  the  inclusion  of  the 
rather  lengthy  computations  required  for  the  radiation  calculation  to  be  carried 
out  without  exceeding  the  high-speed  storage  capacity  of  an  IBM  7090.  The 
synoptic  scale  model  presented  in  this  document  would  have  to  carry  out  such 
computations  through  the  use  of  slower  access  storage. 


T  »’  'TT-t — I'f  ?tf  f  I 


'T'? 


The  computation  of  the  infrared,  radiative -flux  convergence  within  the 
boundary  layer  is  not  independent  of  the  distribution  of  water  vapor,  tempera¬ 
ture,  and  cloudiness  in  the  free  air.  For  short-period,  single-station  forecasts, 
relatively  accurate  estimates  of  these  free  air  distributions  are  obtained  by 
assuming  that  the  initially  observed  state  does  not  vary  appreciably  during  the 
forecast  interval.  This  assumption  would  not  be  very  appropriate  in  a  synoptic 
scale  prediction  scheme.  Moreover,  presently  available  numerical  prediction 
models  for  the  free  air,  in  their  current  form,  do  not  provide  predictions  of  the 
required  parameters. 

In  summary,  we  can  point  to  the  availability  of  procedures  for  computing 
the  influence  of  the  infrared  radiative  flux  which  present  no  formal  obstacle  for 
incorporation  in  our  model.  Oh  the  other  hand,  the  unavailability  of  certain 
required  data  and  the  existence  of  computational  limitations  put  an  adequate  treat¬ 
ment  of  this  process  beyond  our  scope  at  present. 

With  reference  to  the  precipitation  process,  and  more  explicitly,  with  regard 
to  the  evaporation  of  precipitation  falling  into  the  boundary  layer  from  above,  we 
again  are  handicapped  by  the  absence  of  suitable  predictions  from  free  air  models. 
Studies  by  Noguchi  [29],  Caplan  [5]  and  Syono  and  Takeda  [40]  might  be  worked  into 
a  suitable  basis  for  modeling  this  phenomena  if  large-scale  precipitation  rates 
were  provided  by  a  free-air  model.  In  conclusion,  we  point  out  the  potential  for 
improving  our  current  treatment  of  the  physics  of  the  transition  layer.  These 
improvements  may  be  achieved  if  a  suitable  (i.e.  moist)  free-air,  prediction  model 
is  developed. 
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IV 


THE  CONTACT-LAYER  EQUATIONS 

12.  Introduction 

The  transition-layer  equations  require  for  their  solution  the  specification 
of  the  heat  and  moisture  flux  at  the  lower  boundary  of  the  transition  layer.  Also 
required  at  that  level  are  the  horizontal  wind  components  and  the  mixing  coeffi¬ 
cients.  These  conditions  are  provided  in  this  model  by  means  of  a  series  of 
relations  which  are  derived  from  a  number  of  assumptions  regarding  the  structure 
of  the  air  layer  in  close  proximity  to  the  ground.  This  layer  is  referred  to  in  this 
document  as  the  contact  layer.  It  corresponds  to  Lumley  and  Panofsky’s  surface¬ 
boundary  layer  (25],  and  Estoque’s  constant-flux  layer  [13].  There  are  two  cogent 
reasons  for  introducing  this  layer  into  a  synoptic-scale,  boundary-layer  model. 

First,  from  the  computational  viewpoint,  the  accurate  use  of  linear  inter¬ 
polation  in  approximating  derivatives  requires  that  the  parameter  being  differ¬ 
entiated  be  distributed  almost  linearly  between  information  levels.  Near  the 

3 

ground,  atmospheric  variables  tend  to  be  non-linearly  distributed  with  height  . 
Thus,  if  the  prediction  equations  derived  in  Section  III  were  to  be  applied  down 
to  the  ground,  considerations  of  numerical  accuracy  would  require  the  use  of  a 
large  number  of  information  levels  in  this  region,  e.g.,  Fisher  and  Caplan  [15] 
applied  their  prediction  equations  (which  are  similar  to  those  in  Section  III  within 
the  contact  layer,  and  this  required  using  a  vertical  grid-mesh  length  as  small  as 
2  m,  and  concentrating  within  the  lowest  100  m,  twelve  of  the  total  of  twenty-three 
grid  points  used  to  cover  the  full  depth  of  1078  m.  Such  a  disposition  of  informa¬ 
tion  levels  is  disproportionate  for  use  in  the  study  of  synoptic-scale,  boundary- 
layer  processes.  vVhen  the  contact-layer  hypotheses  are  introduced  [13,  30],  one 
may  employ  the  consequent  saving  in  information  levels  to  achieve  greater  refine¬ 
ment  within  the  transition  layer  or  in  a  variety  of  other  ways.  In  our  case,  this 


For  example,  consider  the  validity  of  the  logarithmic  wind  profile  in  a  neutrally 
stratified  boundary  layer. 
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saving  has  been  used  to  study  the  boundary- layer  in  three  space  dimensions. 

The  second  reason  fbr  introducing  the  contact  layer  into  the  model  involves 
physical  considerations.  The  principal  objective  sought  in  application  of  the 
model  derived  here  is  improved  accuracy  in  the  prediction  of  large-scale,  low¬ 
cloudiness.  We  anticipated  the  likelihood  that  turbulent  heat  and  vapor  exchange 
between  the  air  and  ground  would  play  a  part  in  achieving  this  goal.  However, 
it  was  considered  unnecessary  to  achieve  great  accuracy  in  computing  these 
quantities.  Indeed,  if  it  were  true  that  great  accuracy  was  required,  there  would 
be  little  hope  for  success  in  the  application  of  the  model.  One  of  the  most  difficult 
aspects  of  the  micrometeorology  of  the  surface-contact  layer  is  in.  the  specification 
of  the  surface  temperature.  In  a  recent  study,  Elliott  [8]  raises  a  question  as  to 
the  existence  of  a  unique,  ground-air  interface  temperature.  This  question,  and 
others  which  arise  in  a  detailed  study  of  the  interface,  are  of  considerable  importance, 
but  must  be  considered  beyond  the  scope  of  our  present  investigation.  The  contact- 
layer  hypotheses  provide  a  means  for  circumventing  the  requirement  for  a  detailed 
prediction  in  this  region  of  very  complex  physical  interactions,  and  for  obtaining, 
at  least  qualitatively,  correct  estimates  of  the  quantities  needed  to  solve  the 
transition-layer  equations. 

13.  The  Contact- layer  Hypotheses 

By  analyses  such  as  those  given  by  Lumely  and  Panofsky  [25],  one  can 
demonstrate  the  rationale  for  expecting  the  region  near  the  ground  to  be  character¬ 
ized  by  a  small  variation  in  magnitude  of  the  vertical  flux  of  momentum  and  heat. 

To  the  extent  possible  with  instruments  of  limited  accuracy,  micrometeorological 
observations  of  the  eddy  flux  of  heat  and  momentum  support  these  analyses. 

Thus  the  experimental  background  is  present  for  the  working  hypothesis  that,  within 
a  thin  layer  of  air  adjacent  to  the  ground,  the  eddy  flux  may  be,,  fconsidered  invariant 
with  respect  to  height. 

This  assumption  is  a  basis  for  the  development  of  a  similarity  theory  for 
the  structure  of  the  atmosphere  within  this  layer.  Accounts  of  the  derivation  of 
such  theories  are  given  in  Monin  [28],  Priestley  [34],  and  Lumley  and  Panofsky 
[25].  Our  use  of  the  results  of  these  theoretical  investigations  centers  in  the 
derivation  of  formulas  for  evaluating  the  heat  and  vapor  flux  at  the  upper  boundary 
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of  the  contact  layer  ,  We  will  employ  the  Richardson’s  number  as  the  stability 
parameter,  and  neglect  the  difference  between  the  mixing  coefficients  for  heat, 
momentum,  and  vapor. 

The  assumption  that  the  vertical-eddy  fluxes  of  momentum,  heat,  and  vapor 
are  independent  of  height  in  the  contact  layer  permits  one  to  write 


where; 

K  is 
s  is 
T  is 
p  is 
T  is 
y  is 
H  is 

c  is  the  specific  heat  of  air  at  constant  pressure  (0.239  cal  gm  ^  deg  ^), 

P 

q  is  the  specific  humidity, 

-2-1 

Q  is  the  eddy  vapor  flux  (gm  cm  sec  ), 

and  the  quantities  u^,  6^,  and  q^  are  constants  with  the  dimensions  of  velocity, 
temperature,  and  specific  humidity  respectively. 

The  results  of  similarity  theory  and  numerous  empirical  studies  indicate 
that  there  are  two  principal  regimes  which  occur  in  the  contact  layer.  These  are 


^  ^  Z  , 
az  “  P  =  ■ 

K  (—  +  yl  =  —  =  u  9  , 
^dz  c  p  *  * 


(IV-1) 

(IV-2) 


T,  M  “9. 

8z  p  * 


(IV-3) 


the  mixing  coefficient  for  momentum,  heat  and  vapor  (cm^  sec  ^), 

the  horizontal  wind  speed  (cm  sec  ^), 

the  eddy  shear  stress  (gm  cm  ^  sec  ^), 

-3 

the  air  density  (gm  cm  .), 
the  air  temperature  (deg), 

the  adiabatic  lapse  rate  of  temperature  (1°C/100  m), 

-2  -1 

the  eddy  heat  flux  (cal  cm  sec  ), 


-1 


A  unique  upper  boundary  does  not  exist, 
depth  of  the  contact  layer  which  is  a  value 
[25], 


It  is  our  practice  to  choose  50  m  as  the 
suggested  by  Lumley  and  Panofsky 
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denoted  by  the  names;  forced  convection,  and  free  convection.  In  the  forced- 
convective  regime,  the  energy  of  the  turbulent  motions  is  derived  principally 
from  the  kinetic  energy  of  the  mean  flow.  In  the  free-convective  regime,  the 
energy  of  the  turbulence  arises  principally  from  the  bouyancy  forces  associated 
with  the  unstable  density  stratification  produced  in  the  layer  by  surface  heating. 

For  the  forced-convective  regime,  it  follows  [30]  from  the  work  of  Monin 
and  Obukhov  that  the  mixing  coefficient  may  be  expressed  in  terms  of  the  parameters 
of  the  mean  flow  by 

2 

K  =  [kz(l  -  /3R.)]  —  ,  (IV-4) 


where: 


k  is  von  Karman’s  Constant, 

/3  is  a  constant,  and 

R.  is  the  Richardson’s  number. 

1 

The  Richardson’s  number  is  a  dimensionless  ratio  between  the  actions  of  bouyancy 


and  inertia  and  is  defined  by  the  relation, 

R  =  e/H.  J/  [M]' 

i  e  [dz  ' dz 


(IV.-5) 


where  the  ratio  g/ 6  is  the  coefficient  of  thermal  expansion,  g  is  the  acceleration 
of  gravity,  and  6  is  a  mean  temperature  of  the  fluid.  The  term  (1  -  jSR^)  is  the 
analytic  representation  (by  the  first  terms  of  a  power-series  expansion)  of  the 
universal  stability  function  which  is  postulated  to  exist  in  similarity  theory. 

In  the  free-convective  regime,  the  mixing  coefficients  no  longer  depend  upon 
the  stress.  Again,  a  combination  of  similarity  theory  and  numerous  empirical 
investigations  indicates  that,  in  the  free-convective  regime,  the  mixing  coefficients 
may  be  expressed  as  functions  of  the  mean  state  by  the  relation 
o  r  „  ~  n  1/2 


K  =  \z 


2r^  ,aT  _g_ 

F  ‘az  c 


(IV-6) 


The  absolute  value  signs  are  used  because  the  existence  of  the  free-convective 
regime  depends  upon  the  existence  of  a  negative  static  stability.  The  coefficient, 
\,  is  found  to  have  a  value  near  0.9  [34]. 
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Examination  of  observations  (made  under  conditions  with  differing 
Richardson’s  numbers)  of  the  heat  flux  indicates  that  a  relatively  sharp  transi¬ 
tion  between  free-  and  forced-convective  regimes  occurs  about  the  Richardson- 
number  value  of  -.03  [34].  The  importance  of  modeling  such  a  transition  lies  in 
the  different  values  of  the  heat  flux  computed  by  means  of  the  formulas  derived 
for  the  alternate  regimes.  This  point  will  subsequently  be  elaborated. 

If  we  introduce  the  value  for  the  forced-convective,  K  from  Eq.  IV-4,  into 


Eqs.  IV-1,  2,  and  3,  we  may  write 


9z  kz 


1  +  %  - 
e  u. 


az  '  kz  I  e  2 

L  U. 


(IV-7) 

(IV-8) 


dz  kz 


(IV-9) 


Now  consider  the  integration  of  these  equations.  The  wind  speed,  s,  should  vanish 
at  z  =  0  if  the  no-slip  boundary  condition  is  applied.  The  irregularity  of  natural 
terrain  does  not  admit  of  a  unique  level  of  no  slip  and,  in  accordance  with  aero¬ 
dynamic  practice,  a  roughness  length,  z^,  has  been  introduced  by  micromete¬ 
orologists  to  account  for  this  fact  [39].  The  wind  speed  is  then  assumed  to  vanish 
at  that  level.  Integration  of  Eq.  IV-7  then  yields 


s  (z) 


—  §S. 


"  7,  0  -  V’ 


(IV-10) 


where  we  note  that  for  neutral  stratification,  9^  is  zero  and  thus,  Eq.  IV-10  reduces 
to  the  well-known  logarithmic  wind-profile  formula.  The  linear  term  in  this 
expression  will  produce  the  typical  concavity  or  convexity  of  the  wind-speed  profile 
(plotted  on  a  logarithmic  height  scale)  associated  with  non-neutral  thermal 
stratification  [39]. 

As  pointed  out  by  Lumley  and  Panofsky  [26],  the  constant- flux  hypothesis 
as  expressed  in  Eq.  IV-8  is  of  doubtful  validity  below  z  =  1  m  due  to  the  relatively 
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strong  convergence  of  the  radiative  flux  below  that  height.  Consequently,  denoting 

the  level  of  the  instrument' shelter  as  z  =  z.  ^  1  m,  we  integrate  Eq.  (IV-8)  using 

1 

the  boundary  condition  T  =  T.  at  z  =  z..  The  result  is 


I  z  1  j3g  i  * 

T(Z)  =  T.  -  v(z  -  z.)  .  ITl-'l-h  T'-|  -  V 


(IV-11) 


For  the  sake  of  compatibility  with  Eq.  IV-11,  we  will  specify  the  boundary  condition 
for  Eq.  iy-9  at  z  =  z^  also.  The  integration  yields 


k  z.  0 


q(z)  =  q(Zi)  +  T7  In'  —  +  ^  (z  -  z  ) 


(IV-12) 


An  analogous  set  of  profile  equations  may  be  derived  from  the  equations  of 
constant  flux  by  use  of  the  free-convective-regime  formula  for  the  mixing  coeffi¬ 
cient  (Eq.  IV-6).  The  boundary  conditions  are  applied  at  the  same  levels  as  in  the 
case  of  forced  convection.  One  finds  that  the  profile  equations  may  be  written  as 

3“.^  /„l-l/3r-l/3  -l/sl 


-1/3  -1/3 

T(z)  =  T  -  y(z  -  z  )  -  — - —  I  [z  -  z  ], 


3u,q,  1/3 

q(z)  =  q.  -  77^  w  [z  -  z  ]. 


(IV-14) 

(IV-15) 


14.  The  Empirical  Formulas 

Estoque  and  Pandolfo,  in  their  previously  cited  work  with  boundary-layer 
prediction  models,  have  employed  the  contact-layer  profile  relations  derived  above. 
Their  use  of  them  differs  in  an  important  respect  from  ours.  The  difference  is 
necessitated  by  our  omission  of  a  prediction  equation  for  the  horizontal  wind  in 
the  transition  layer.  This  omission  results  from  our  conviction  that  it  would  have 
been  inadvisable  to.  embark  on  this  investigation  with  a  dynamically  active  model. 


especially  in  view  of  the  well-known  inadequacy  of  current  synoptic  wind  observa- 
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tions  for  specifying  the  geostrophic  deviation  . 

The  use  of  the  profile  formulas  in  our  work  will  be  spelled  out  in  the  next 
two  sections.  In  this  section  we  will  discuss  the  empirical  formulas  for  the 
stress  and  geostrophic-deviation  angle  which  form  the  basis  upon  which  the 
subsequent  formulation  is  constructed. 

The  empirical  evidence  which  is  employed  here  is  taken  from  a  paper 
presented  by  Prof.  A.  K.  Blackadar  [3]  before  a  national  meeting  of  the  American 
Meteorological  Society.  This  paper  deals  with  the  steady-wind  distribution  in  a 
neutral,  baroclinic  boundary  layer.  In  presenting  his  results,  Blackadar  used 
observational  data  on  the  relations  between  the  friction  velocity  u^,  the  geostrophic- 
wind  speed,  G,  the  surface  Rossby  number  (G/fz^)  and  the  angle  of  deviation 
between  the  surface  wind  and  the  surface  geostrophic  wind.  The  sources  of  these 
observational  data  are  identified  in  his  paper. 

We  found  that  the  following  expressions  could  be  used  to  approximate  the 
principal  characteristics  of  these  relations: 

u^  =  G  [0.07625  -  0.p0625r),  (IV-16) 

2 

ip  =  0.625  r  -  12.750  r  +  80.625,  (IV-17) 

where 

r  =  logio(^),  (IV-18) 

u^  is  the  friction  velocity,  G  is  the  surface  geostrophic  wind  speed,  f  is  the 
Coriolis  parameter,  and  the  surface  roughness.  The  quantity  tp  is  the  angle 
(in  degrees)  between  the  surface  wind  and  the  surface  geostrophic  wind.  Blackadar ’s 
theoretical  results  indicate  that  the  wind  direction  does  not  vary  appreciably  through 
the  lowest  50  m,  and  this  is  consistent  with  our  use  of  the  constant-flux  hypothesis. 

These  empirical  formulas  are  not  dependent  upon  the  static  stability  of  the 
air  in  the  contact  layer.  There  is  some  theoretical  evidence  quoted  in  Haltiner  and 


5 

In  this  connection,  simply  recall  the  fact  that  the  majority  of  wind  reports  are 
based  on  pilot  balloon  runs  followed  visually. 
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Martin  [19]  that  the  geostrophic-deviation  angle  should  tend  to  be  greater 
in  inversion  conditions  and  smaller  in  lapse  conditions.  According  to 
Blackadar’s  theoretical  results,  the  baroclinicity  of  the  boundary  layer  also 
accounts  for  a  tendency  toward  larger  deviation  angles  when  the  thermal  wind 
direction  is  opposite  to  the  surface  geostrophic  wind  direction  and  vice  versa. 
Thus,  we  note  that  Eq.  IV-17  does  not  take  into  account  two  of  three  relevant 
physical  factors.  However,  there  seemed  to  be  little  point  in  trying  to  model 
these  refinements  until  some  basic  evaluations  of  the  model  are  carried  out. 
This  is  especially  true  because  the  geostrophic  wind  is  not  predicted  by  this 
model;  it  must  be  obtained  from  an  independent,  free-air  prediction  model. 

The  friction-velocity  computation  based  on  Eq.  IV-16  may  also  require 
modification  to  take  into  account  the  existence  of  thermal  enhancement 
(suppression)  of  the  turbulent  exchange  in  lapse  (inversion)  conditions.  In 
some  initial  idealized  experiments,  analogous  to  those  quoted  earlier,  we  found 
that  use  of  Eq.  IV-16  tended  to  yield  as  large  a  negative  heat  flux  values  as 
positive.  This  seemed  sufficiently  unrealistic  to  require  modification  of 
Eq.  IV-16  even  in  our  basic  model.  We  therefore  replaced  Eq.  IV-16  with  the 
following  pair  of  equations; 


u*  =  1.2  G  [0.07625  -  0.00625r],  T"  ^  0 

9z 

u,^  =  0.8  G  [0.07625  -  0.00625  r],  ~  >  0 

oz 


(IV-19) 


Repetition  of  the  idealized  experiments  supported  the  use  of  these  formulas. 
In  practical  application  of  a  model  similar  to  this  one,  we  will  have  to  recognize 
again  the  importance  of  having  available  an  accurate  prediction  of  the  geostrophic 
wind. 


15.  Eddy- flux  Evaluation 

The  technique  employed  in  the  computation  of  the  eddy  fluxes  of  heat  and 
vapor  at  the  interface  (z  =  h)  of  the  transition  and  contact  layers  will  be  explained 
in  this  section.  The  transition-layer  prediction  equations  are  applied  at  the  level 
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z  =  h  to  obtain  values  of  temperature  and  specific  humidity  at  this  level  . 

These  quantities,  together  with  the  time-dependent  boundary  conditions  at 
z  -  z.,  are  used  in  the  profile  formulas  to  transform  them  to  equations  in  the 
quantities  9^  and  q^.  The  empirical  calculation  of  u^^  permits  us  to  consider  it 
as  a  known  quantity  in  the  profile  equations  for  temperature  and  vapor. 

The  equations  for  the  fluxes  will  be  evaluated  for  the  forced-convective 
regime.  We  will  show  that  these  formulas  fail  when  the  value  of  Richardson’s 
number  is  nearly  equal  to  the  critical  Richardson’s  number  discussed  by 
Priestly  [34].  A  transition  to  the  free-convective  formulas  is  then  indicated, 
and,  therefore,  the  flux  formulas  for  free  convection  are  derived.  We  conclude 
this  section  with  a  summary  of  the  computational  formulas. 

We  denote  the  values  of  temperature  and  specific  humidity  predicted  at 
the  level  z  =  h  by  the  subscripted  quantities,  T^^  and  q^^,  respectively.  Consider 
first  Eq,  IV-11,  the  equation  for  the  temperature  profile  in  forced  convection. 

If  it  is  evaluated  at  z  =  h,  the  equation  becomes  a  quadratic  in  6^.  The  two  roots 
of  this  equation  are  readily  obtained.  Upon  inspection,  and  realizing  that  6^ 
must  vanish  when  the  lapse  rate  is  neutral,  the  appropriate  root,  or  solution,  is 
found  to  be 


^  -0  (h/Z.)  j 

2/3  ‘  g  ■  k(h  -  Z.)  \ 


(IV- 20) 


M  £ 

2  ■  F 
u* 


T.  -  T  \  ;  In  h/Z.  \  -2 

_i _ h  _£_1(  _ L_ 

h  -  Z.  C  I  1  k(h  -  Z.) 

1  P  1 


Consideration  of  Eq.  IV-11  evaluated  at  z  =  h  shows  that  it  is  readily  evaluated 


for  q^.  The  result  is 


q*  =  (Qh  -  q^) 


Inh/Z. 


-2  |  -  Z 


1  -1 
i)  - 


{IV- 21) 


We  do  not  believe  that  any  real  problem  arises  from  the  dual  character  of  the  level 
z  =  h  in  our  use  of  both  contact-  and  transition-layer  formulas  at  this  level. 

Estoque  [13]  uses  a  device  to  circumvent  this  question,  but  the  quantitative  cal¬ 
culations  are  sensibly  unmodified. 
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The  heat  flux  and  vapor  flux  within  the  contact  layer  in  forced  convection  are 
then  obtained  straight-for’ivardly  by  multiplying  6^  and  (computed  by  the 
formulas  above)  by  u^^  (computed  from  the  empirical  formulas).  We  then  note 
that  Eq.  IV-20  yields  a  complex  value  of  0^  whenever  the  temperature  gradient 
satisfies  the  following  inequality; 


C 

P 


—  i 
^  4/3  g 


In  h/Z.  I  2 
k(h  -  Z.) 


(IV-22) 


Thus,  Eq.  IV-20  can  only  be  applied  with  a  temperature  gradient  at  least  as 
stable  as  the  limiting  value  given  by 


g  4/3 


(IV -23) 


At  this  limiting-temperature  distribution,  the  heat  flux,  H,  computed  from  the 


forced-convective  formula  is 


ec  pu^  In  h/Z. 
H  =  -  pCp  u^e^  =  -  Z.) 


(IV- 24) 


In  the  forced-convective  regime,  the  Richardson’s  number  can  be  written 


/3gkz  0*  >  -1 


R.  =  -  1  -  1  + - 

1  /3  L  1 


(IV-25) 


Now,  at  the  limiting-temperature  gradient  given  in  Eq.  IV-23,  the  Richardson’s 


nuniber  may  be  evaluated  using  Eq.  IV-24.  The  result  is 
,  r  In  h/Z,  .-It 


R.  =  -  1  -  1  -  - 

1  /3  2 


-M  i’ 


(IV- 26) 


This  Richardson’s  number  may  be  evaluated  for  various  heights,  z,  and  values 
of  the  parameter  /3  (See  Table  I  )..  According  t,o  Priestley  [34],  the  transition 
from  forced-to  free-convection  is  associated  with  a  Richardson  number  of  -0.03 
or  -0.02,  rpeasured  at  ;  z  =  1.5  m.  Pandolfo  [30]  (also  see  Estoque  [13] 
recommends  the  use  of  /3  =  3.0  in  the  forced-convective  formulas  on  the  basis 
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of  his  analysis  of  profiles  in  inversion  conditions..  From  Table  I  we  see  that 

I 

at  z  -  l.'s  m  with  /3  =  3,  the  Richardson’s  number  associated  with  the  limiting 
temperature  distribution,  Eq.  lV-23,  agrees  with  Priestley’s  value. 

TABLE  I 

RICHARDSON’S  NUMBER 

(associated  with  limiting  temperature  distribution  between  1  and 
50  meters  for  forced-convective  heat  flux  equation)  as  a 
function  of  the  parameter,  j3,  and  the  height  z(m)  at  which  R.  is 
measured  in  a  layer  of  constant- flux  of  total  depth  50  m. 


i3  X 

"■>  .5 

1.0 

1.5 

2.0 

2.5 

3.0 

1 

-.021 

’-.039 

-.063 

-.087 

-.111 

-.135 

2 

-.010 

-.019 

-.032 

-.044 

-.056 

-.007 

3 

-.007 

-.013 

-.021 

-.029 

-.037 

-.045 

4 

-.005 

-.009 

-.016 

-.022 

-.028 

-.034 

5 

-.004 

-.008 

-.012 

-.017 

-.022 

-.027 

10 

-.002 

-.004 

-.006 

-.009 

-.011 

-.014 

Because  the  forced-convective  heat-flux  formula  fails  to  supply  a  useful 
result  on  the  lapse  side  of  the  limiting  temperature  distribution,  one  may  there 
apply  the  free  convective  formulas.  It  is  of  some  interest  to  examine  the  pos¬ 
sibility  of  enforcing  continuity  in  the  heat  flux  computed  from  the  forced-and 
free-convective  equations  at  the  limiting  temperature  distribution.  In  the  first 
instance,  the  profile  formula,  Eq.  IV-14,  may  be  evaluated  at  z  =  h  and  solved 


Now,  at  the  limiting  temperature  distribution  given  by  Eq.  IV-23,  the  free- 


convective  heat  flux  is  given  by  Eq.  IV-27  as 
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_l/3  /31-3/2  (Inh/Zp 

3^/"  ‘  ■  J  l4^k^h■-  Z,)_ 


2  -13/2 


(IV-28) 


For  continuity  in  the  computed  heat  flux  at  the  limiting  temperature  distribution, 
Eqs.  IV-28  and  24  should  be  equal.  This  equality  then  implies  a  relationship 
between  )3  and  which  after  some  manipulation  may  be  solved  for  \  to  give 


,  3/2  -1/3  , -1/33/2  ,  ^  1/2 

4(3)  '  (Z  -  h  '  )  '  (h  -  Z  )  '  / 

\  =  k'  - ^ - z - ^ 

(Inh/Z.) 

1 

Again,  using  =  1  m,  h  =  50  m  and  k  =  0.38,  we  find  that 


(IV-29) 


\  =  0.85 


(IV-30) 


Priestley  [34]  quotes  0.9  as  the  most  likely  value  of  \,  whereas,  according  to 
Eq.  IV-30,  we  note  that  for  )3  =  1,  \  =  .85;  for  ^  =  2,  X  ~  1.2;  and  for  j3  =  3, 

\  =  1.5.  The  values  )3  =  2  and  \  =  1.2,  although  not  the  precise  values  indicated 
by  the  experimental  evidence,  are  considered  to  be  sufficiently  accurate  for  our 
work.  If  the  heat  flux  in  free  convection  is  computed  from  Eq.  IV-27,  it  is  simple 
to  evaluate  the  vapor  flux  from  Eq.  IV-15.  We  will  next  recapitulate  the  procedure 
for  computing  the  heat  and  vapor  fluxes  in  the  contact  layer. 

The  first  step  is  to  compute  u,^,  from  the  appropriate  formula  of  the  pair 
denoted  as  Eq.  IV-19.  One  may  then  determine  if  the  free-  or  forced-convective 
formulas  are  appropriate  by  evaluating  the  left-  and  right-hand  sides  of  Eq.  IV-23. 

If  forced  convection  is  indicated,  one  solves  Eq.  IV-20  for  6^,  and  then 
solves  Eq.  IV-21  for  q,|^.  The  forced-convective  fluxes  are  then  computed  from 
u„,  and  q„. 

If  free  convection  is  indicated,  the  heat  flux  is  computed  from  Eq.  IV-27. 
Knowing  Eq.  IV-15  may  be  evaluated  for  the  vapor  flux. 


16.  Formulas  for. the  Horizontal  Wind  Components  and  Mixing  Coefficients 
In  Subsection  13  we  presented  the  following  expressions  for  the  mixing 
coefficient  in  forced  and  free  convection,  respectively. 


K  =  [kz(l  -  )]^  |J, 

1  dZ 


(IV-31) 
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K  =  K(h)  (IV-38) 

The  boundary  values,  U  and  V,  required  for  computing  the  horizontal  wind 
in  the  transition  layer  are  derived  as  follows. 

In  forced  convection,  Eq,  IV-10  yields 

s(h)  =  —  In  ^  (h  -  Z  )  (IV-39) 

k  Zq  u^  0 
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For  free  convection  Eq.  IV-13  yields 


3u. 


s  (h)  =  - 


2/3  1/3 


>-1/3 

£ 

-  z 

0  1 

0 

(IV-40) 


Let  the  geostrophic  wind  components  at  z  =  h  (evaluated  from  Eqs.  III-4a 


and  b)  be  denoted  by  U  and  V  .  One  may  then  fix  the  angle 

S  S 


0  =  tan 

g 


-1 


V 

L^g  J 


(IV-41) 


in  degrees.  Combining  0  and  the  deviation  angle,  tp  computed  from  Eq.  IV-17, 

g 

we  have 


0  =  0g  +  ^, 


(IV-42) 


and 


U  =  s  (h)  cos  0 
V  =  s  (h)  sin  0 . 


(IV-43) 


17.  Contact-layer  Boundary  Conditions 

This  model  has  been  designed  so  that  for  the  transition-layer  prediction 
equations  the  lower-boundary  conditions  are  derived  from  formulas  obtained 
from  the  equations  valid  in  the  contact  layer.  The  latter  set  of  equations  also 
requires  the  availability  of  boundary  conditions  for  their  evaluation.  One  set  of 
these  conditions  is  supplied  by  requiring  continuity  in  the  dependent  variables  at 
the  interface  between  the  two  layers.  There  remains  the  need  to  provide  a 
means  for  computing  the  temperature  and  specific  humidity  at  the  level  of  the 
instrument  shelter,  z  = 

We  will  consider  first  the  procedure  for  obtaining  the  temperature  at  the 
instrument  shelter.  The  primary  cause  for  variations  in  this  temperature  is  the 
variation  in  the  solar  radiation  received  at  the  ground.  Secondary  factors  involve 
terrestrial  and  atmospheric  radiation.  All  of  these  factors  are  influenced  by  the 
constitution  bf  the  atmosphere  (clouds,  water  vapor,  carbon  dioxide,  etc.)  and  by 
the  ground  characteristics  (albedo,  thermal  conductivity,  etc.).  Because  of  this 
complex  interaction,  it  did  not  seem  realistic  to  attempt  in  this  model  to  specify 


33 


deterministically  the  temperature  at  the  level  of  the  instrument  shelter.  One 
might  achieve  adequate  accuracy  through  the  use  of  a  combination  of  statistical 
prediction  methods  and  physical  climatology.  This  is  one  of  the  problems  we 
plan  to  study  during  the  coming  year.  In  the  experiments  reported  in  this  paper 
(see  Section  VII),  we  neglected  the  change  in  temperature  due  to  radiation. 

In  addition  to  radiation  processes,  the  temperature  at  the  level  of  the 
instrument  shelter  may  be  expected  to  change  in  response  to  the  advection 
process.  Because  the  instrument- she  Iter  level  is  within  the  contact  layer,  we 
decided  to  omit  advective-temperature  changes  from  consideration  unless  they 
were  greater  than  the  temperature  change  which  would  be  produced  by  a  con¬ 
vergence  within  the  contact  layer  of  some  20  percent  of  the  heat  flux  passing 
through  the  layer.  The  procedure  used  in  the  experiments  (see  Section  yil) 
estimated  the  temperature  change  due  to  advection  at  the  top  of  the  contact  layer 
and  applied  this  change  to  the  instrument  shelter  whenever  it  exceeded  the 
criterion  given  above. 

In  determining  the  specific  humidity  at  the  level  z  =  z.,  we  assumed  the 
relative  humidity  was  specified  a  priori.  This  value,  coupled  with  the  scheme 
given  above  for  obtaining  the  temperature  and  an  assumption  that  the  air  pressure 
was  known  at  z  =  z^,  is  sufficient  for  determining  the  specific  humidity  at  z  =  z^. 
Here  again  we  believe  that  statistical  prediction  methods  utilizing  those  pre¬ 
dictors  available  within  the  model  could  provide  an  improved  technique  for 
prescribing  the  relative  humidity. 
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V 


COMPUTATIONAL  FORMULATION  OF  MODEL 


In  the  preceding  sections,  the  various  equations  which  govern  the  behavior 
of  our  model  of  the  atmospheric  planetary-boundary  layer  were  formulated.  We 
now  describe  the  computational  scheme  by  which  the  model  equations  are  to  be 
applied  to  the  prediction  problem.  The  differential  equations,  region  of  integra¬ 
tion,  and  time  are  made  discrete  by  the  introduction  of  a  grid  mesh  of  points. 

The  region  of  integration  will  be  of  limited  horizontal  extent..  For  sim¬ 
plicity,  we  will  assume  that  the  grid  points  in  the  (x,  y)-plane  are  equally  spaced 
a  distance  d  apart.  The  grid  points  along  each  vertical  will  be  spaced  at  unequal 
distances,  (Az).,  apart.  The  density  of  grid  points  in  the  vertical  will  diminish 
with  increasing  distance  above  the  ground.  The  time  coordinate  will  be  divided 
into  increments  of  equal  duration.  At.  We  will  assume  that  the  region  may  be 
covered  by  (L  +  1)  grid  points  in  the  x-coordinate,  (M  +  1)  grid  points  in  the 
y-coordinate,  and  K  grid  points  in  the  z-coordinate.  We  will  use  the  following 
notation  to  represent  a  function’s  value  at  a  grid  point; 


where 


.l.n 

m,  k 


=  f(x=x^,  y=y_,  z=z,^,  t=t^) , 


m 


x^  =  (1-1)  d, 
Ym  =  («i-l)d, 


k 

i=l 


tn  =  T  +  {n-l)At, 


n 

1  =  1,  ...,  L+1 ; 
m  =  1,  ...,  M+1 ; 

k  =  1,  ...,  K; 
n  =  1,  ...,  N+1  ; 


(V-1) 


and  T  is  the  time  of  the  initial  data  measured  from  some  suitable  reference  time. 

If  we  wish  to  refer  to  the  value  of  a  function  at  some  particular,  but  un¬ 
specified  value  of  X,  y,  z,  or  t,  we  will  use  a  Greek  letter  in  place  of  a  Roman 

symbol.  Thus,  the  value  of  a  function  at  a  particular  time  step  may  be  denoted 
.1.^^  ■ 


by  f 


m,  k* 
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18.  The  Difference  Equations 


The  prediction  equations  for  the  transition  layer  will  be  put  into  a  partially- 
implicit  difference  form.  This  partially-implicit  scheme  has  the  major  advantage 
of  being  computationally  stable  for  relatively  large  time-step  increments.  One 
must  pay  for  this  advantage  by  solving  the  set  of  simultaneous  linear  equations 
which  result  from  the  difference  equations.  A  simple  computer-oriented  solution 
of  these  algebraic  equations  may  be  obtained  by  the  method  of  Gaussian  Elimina-  ■ 
tion.  Before  presenting  the  difference  equations,  we  discuss  the  form  taken  by 
the  various  differential  terms. 

The  time  derivative  in  all  the  equations  is  approximated  by  a  forward 
difference  over  an  interval,  At.  We  write 


,l,n+l  l,n 


at 


m,  k 


m,  k 


(V-2) 


At 


It  is  to  be  noted  that  computational  stability  is  insured  by  the  partially-implicit 
difference  scheme  for  relatively  large  values  of  (At)  (see  subsection  19),  but 
that  in  order  to  insure  a  good  approximate  solution,  the  value  chosen  for  At  must 
not  be  too  large. 

The  horizontal  advection  terms  are  approximated  by  the  upwind  technique 
used  by  Fisher  [14]  and  discussed  in  Forsythe  and  Wasow  [16],  This  procedure 
is  outlined  in  Appendix  II  for  the  various  advection  quantities.  The  one-sided 
derivative  approximation  is  theoretically  of  lower-order  accuracy  than  that 
obtainable  with  a  centered  difference  approximation,  but  from  a  practical  view¬ 
point  the  method  is  in  better  accord  with  the  concept  of  advection  as  a  transport 
process  dependent  only  upon  upstream  values  of  the  transported  quantity,  and 
will  frequently  involve  smaller  truncation  errors. 

The  diffusion  terms  and  the  vertical  advection  terms  are  approximated  by 
centered  space  derivatives.  The  mixing  coefficients  and  vertical  velocities  are 
evaluated  at  the  current  time  step,  whereas  the  quantities  being  diffused  and 
advected  are  evaluated  at  the  subsequent  time  step.  This  procedure  introduces 
the  so-called  implicit  character  into  the  difference  equations.  The  vertical 
coordinate  index  on  the  mixing  coefficients  is  that  of  the  grid  point  at  the  top  of 
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the  interval  for  which  the  coefficient  is  evaluated. 

The  difference  equations  approximating  the  prediction  equations  are  written  out 
in  Appendices  III,  IV,  and  V.  In  order  to  demonstrate  the  method  of  solution'of 
the  difference  equations,  it  is  adequate  to  treat  just  one  equation;  the  heat  trans¬ 
fer  equation  has  been  selected  for  this  purpose.  Using  the  equations  outlined  in 
Appendix  III,  cne  may  rewrite  these  equations  in  the  form 


a"-:  .  b'’,"  .  c"',’'  (k=l . K-D- 

p,k  p,k  ^l,k  M,k+1  p,k 


(V-3) 


The  various  coefficients  in  Eq.  V-3  are  spelled  out  in  Appendix  VI.. 


The  system  of  (K-1)  algebraic  equations  represented  by  Eq.  V-3  may  be 
solved  for  the  (K-1)  values  of  (k=l,  ...,  K-1)  by  the  method  of  Gaussian 

Elimination  which  does  not  involve  any  iteration.  The  various  coefficients  in 
Eq.  V-3  are  all  computed  from  quantities  known  at  t  =  t^.  The  upper -boundary 
temperature  ^  is  known  at  t  =  t^,  because  it  is  a  prescribed  boundary 

condition.  The  other  prediction  equations  are  susceptible  to  an  identical  solution 
procedure. 


19;  Analysis  of  Computational  Stability 

The  choice  of  difference  equations,  indicated  in  Appendices  III,  IV,  and  V, 
was  made  to  permit  the  use  of  a  relatively  long  time  step.  It  is  well  known  that 
the  solution  of  a  difference  equation  may  exhibit  an  improper  exponential  growth 
if  the  ratio  of  the  time  and  space  intervals  does  not  satisfy  certain  criteria.  These 
computational  stability  criteria  are  derived  from  consideration  of  linearized  ver¬ 
sions  of  the  difference  equations  in  which  the  coefficients  are  assumed  to  be 
constants.  The  intuitive  justification  for  regarding  such  criteria  as  applicable  to 
the  more  general  difference  equations  is  presented  in  Section,  IV  of  the  paper , 
by  Richtmyer  [36]. 

Because  the  various  equations  employed  in  this  model  are  of  similar  form, 
we  will  derive  the  computational  stability  criteria  for  a  generic  case.  For 
simplicity,  only  the  principal  steps  in  the  analysis  are  presented  here,  and  we 
retain  only  one  of  the  horizontal  coordinates. 
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,  iKAz  -iKAz  ,  .  ,  i/cAz  -iKAz  „ ,  ^  1  ,  (V-9) 

l  +  p[e  -e  ]  +  y[e  +e  -2]  ' 

which  is  a  necessary  condition  for  the  difference  equation  to  possess  bounded 

7 

solutions  of  the  form  of  Eq.  V-6  for  n  . 

2  2  2  2 

a((7  -  1)  (1  -  cos  yAx)  2y(l  -  cos /cAz)  +  2y  (1  -  cos  kAz)  +  2p  sin  kAz 

(V-10) 

Inspection  shows  that  the  right  hand  side  of  Eq.  V-10  is  non-negative, 
whereas  the  left  hand  side  will  be  non-positive  if  ct  1.  Thus,  we  conclude  that 
the  difference  equation  will  have  computationally-stable  solutions  provided  that. 

Ay 

At  ^  {-^).  ,  (V-11) 

This  requires  that  the  time  step  not  exceed  the  space-grid  mesh  divided  by 
the  largest  horizontal  velocity  which  is  likely  to  occur.  If  we  assume  that  Ax  is 
100  km  and  U  is  40  m  sec  then  At  may  be  as  large  as  40  minutes. 

A  theorem,  due  to  Lax,  which  states  the  equivalence  of  stability  and  con¬ 
vergence  for  properly  posed  initial-value  difference  equations  is  proven  in 
Richtmyer’s  book  [36].  This  implies  for  our  problem  that,  if  At  is  taken  suffi¬ 
ciently  small  while  the  condition  Eq.  V-11  is  maintained,  the  solution  obtained 
from  the  difference  equation  will  be  close  to  the  solution  of  the  differential  equation. 
This  difference  between  the  solutions  of  the  difference  and  differential  equations 
can  be  ascribed  to  truncation  error.  It  would  seem  necessary  then,  to  experiment 
with  various  choices  of  At  which  satisfy  Eq.  V-11  to  arrive  at  an  optimum  value 
relative  to  the  two  desirable  features: 

(a)  Minimal  computation  time  for  a  forecast  of  fixed  duration,- 

(b)  Minimal  truncation  error  relative  to  the  observational  accuracy  of  the 
data  employed. 

7 

n  —  implies  either  (a)  for  At  fixed,  the  duration  of  the  prediction,  T,  — 
or(b)  for  T  fixed.  At  approaching  zero. 


Such  experiments  may  be  carried  out  in  conjunction  with  investigation  of  the 
relationships  among  stability,  convergence  and/ or  truncation  error  for  finite 
values  of  At,  and  the  space  increments  Az  and  Ax. 

20.  Specification  of  Boundary  Conditions 

In  Section  III,  the  boundary  conditions  necessary  for  the  solution  of  the 
differential  equations  were  enumerated.  These  conditions  apply  equally  to  the 
difference  equations  and  are  of  two  kinds.  At  the  upper  boundary,  the  values  of 
the  temperature,  specific  humidity,  specific  moisture,  geostrophic  wind  com¬ 
ponents,  and  thermal  wind  components  are  required.  At  the  lower  boundary,  the 
fluxes  of  heat  and  vapor,  together  with  the  components  of  the  horizontal  wind, 
are  required. 

The  first  kind  of  data  is  to  be  provided  independently  of  the  boundary  layer 
model.  It  is  our  present  plan  to  specify  these  values  in  tests  of  the  boundary- 
layer  model  by  use  of  the  predictions  made  with  an  Air  Weather  Service  multi¬ 
level  model  or,  when  such  are  not  available,  from  observed  data.  More  particularly, 
the  predicted  values  of  temperature,  specific  humidity,  and  geopotential  height  for 
850,  700,  and  500  mb  at  six-hour  intervals  are  used  with  linear  interpolation  to 
provide  the  needed  quantities  at  2000  meters  above  the  ground  at  each  time  step. 

In  order  to  provide  the  lower -boundary  conditions,  we  make  use  of  the 
various  formulas  outlined  in  Section  IV.  These  require  for  their  evaluation  the 
predictions  made  by  the  boundary-layer  model  at  z  =  h,  as  well  as  the  surface 
temperature  and  relative  humidity  which  are  currently  provided  by  the  schemes 
discussed  in  subsection  17. 
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VI 

EXPERIMENTAL  BACKGROUND 

21.  Preliminary  Experiments 

The  initial  experiments  conducted  with  the  model  equations  were  simplified 
by  using  only  two  space  dimensions,  (one  vertical  and  one  horizontal)  and  the  use 
of  idealized  initial  and  boundary  conditions.  The  variety  of  results  obtained  in 
this  manner  is  not  presented  here,  but  these  results  were  used  to  modify  the 
various  equations  of  the  model.  Some  of  these  instances  of  model  modification 
were  pointed  out  in  the  preceding  text.  In  this  sense,  the  results  presented  in 
Section  VII  of  this  report  have  also  led  us  to  modify  certain  parts  of  the  model. 
The  section  dealing  with  the  use  of  the  free- convective  formulas  is  an  example 
of  such  a  change.  Unfortunately,  we  did  not  have  time  to  re-program  this' 
change  into  the  model  and,  consequently,  the  experiments  reported  in  Section  VII 
do  not  permit  us  to  evaluate  its  significance.  This  state-of-flux  in  the  model 
formulation  indicates  that  the  model  should  be  considered  a  base  model,  similar 
to  a  electronic  engineer’s  breadboard  design.  In  the  following  subsections,  we 
will  discuss  the  subsidiary  analysis  techniques  used  to  process  the  observed  data 
and  to  specify  the  exterior  factors  involved  in  the  model. 

22.  The  Analysis  of  Initial- Temperature  and  Humidity  Data 

We  collected  the  synoptic  upper-air  and  surface  data  (received  via  telet5q3e) 
at  the  weather  station  of  the  TRC  Service  Corp.  during  the  five-day  period, 

6  February— 10  February,  1964.  In  order  to  process  the  radiosonde  observations, 
we  first  used  linear  interpolation  between  the  significant-level  reports  to  obtain 
the  temperature  and  dew  points  at  the  levels  required  in  the  model.  This  step 
was  accomplished  through  the  use  of  the  computer  program,  “Sounding  Construc¬ 
tion,”  written  at  the  UAC  Computation  Laboratory. 

These  data  were  available  at  the  various  radiosonde  observing  stations 
located  within  the  geographical  region  to  which  the  model  was  to  be  applied.  In 
order  to  start  the  numerical  computation  of  the  forecast,  these  data  must  be 
available  on  a  regular  array  of  grid  points.  With  the  intention  of  automating,  and 
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making  objective,  the  analysis  of  the  initial  data,  we  had  programmed  a  routine 
for  interpolating  from  the  data  at  the  radiosonde  stations  onto  the  horizontal 
grid  point  array.  The  Successive  Approximation  Technique  (SAT)  system  [41] 
for  performing  such  an  analysis  was  used.  This  technique  requires  the  pre¬ 
liminary  specification  of  a  first  guess  to  a  grid-point  analysis.  We  attempted 

8 

to  use  this  method  with  two  types  of  initial  guess  .  The  results  were  not  satis¬ 
factory,  principally  due  to  the  absence  of  a  means  of  enforcing  vertical  consistency 
in  the  resulting  grid-point  data  analysis.  We  therefore  decided  to  abandon,  for 
the  time  being,  the  attempt  to  completely  automate  the  analysis. 

The  method  finally  adopted  involved  subjective  analysis.  The  difference  in 
temperature  between  successive  levels  in  the  vertical  was  computed  for  each 
reporting  station.  These  numbers  were  then  plotted  on  a  base  map.  On  each 
map,  isopleths  were  drawn  for  the  values  of  the  vertical-temperature  difference 
(proportional  to  the  static  stability  of  the  layer).  These  isopleths  were  then  used 
to  interpolate  the  temperature  differences  onto  the  grid  points.  The  data  at  the 
grid  points  on  each  map  were  then  used  in  conjunction  with  a  surface-temperature 
analysis  to  reconstruet  the  temperature,  at  each  level,  for  every  grid  point. 

In  analyzing  the  humidity  field,  we  used  the  results  of  the  “Sounding 
Construction”  program  to  compute  the  dew-point  depression  at  each  level. 

These  were  then  analyzed  in  the  manner  used  for  the  temperature  difference. 

The  final  dew  point  distribution  resulted  from  a  combination  of  the  independent 
analyses  of  temperature  and  dew-point  depression. 

23.  The  Analysis  of  the  Geostrophic  Wind 

In  Subsection  5,  we  noted  that  the  geostrophic  wind  was  assumed  to  vary 
linearly  with  the  height  through  the  boundary  layer.  It  is  rather  straight  forward 
to  compute  the  geostrophic  wind  at  the  surface  under  the  assumption  that  it  is 
equivalent  to  that  computed  from  sea-level  pressures.  Similarly,  the  computation 


These  methods  are  given  in  the  monthly  progress  report  for  July,  1964. 
(External  Monthly  Progress  Report,  Contract  No.  11910,  The  Travelers  Research 
Center,  Inc.) 
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of  the  geostrophic  wind  in  the  upper-air  surfaces  of  constant  pressure  may  be 
computed  directly.  In  practice,  we  interpolated  the  surface  pressure  and  upper- 
air  geopotential  anal3'^ses  into  grid-point  values,  and  evaluated,  by  finite  differ¬ 
ence  approximation,  the  geostrophic-wind  components  at  mean  sea  level,  850  mb 

0  0 

and  700  mb.  The  sea-level  components  were  used  to  evaluate  u  and  v  in 

S  ^  H 

Eqs.  III-4a  and  b.  The  determination  of  the  geostrophic-wind  components,  u 

H  ^ 

and  V  ,  valid  at  the  upper  boundary  of  the  model-boundary  layer,  was  done  by 

S 

linear  interpolation  between  the  850-  and  700-mb  geostrophic-wind  components. 
This  interpolation  was  based  on  the  reported  height  of  these  pressure  levels 
above  the  station. 

24.  Terrain  Height  and  Surface  Roughness 

Figure  1  presents  the  distribution  of  the  height  above  mean  sea  level  of 
the  grid  points  used  in  the  computation  discussed  in  Section  VII.  The  heights 
are  plotted  in  meters,  with  isopleths  drawn  at  100  meter  intervals.  The  data 
used  in  arriving  at  these  figures  were  taken  from  Berkofsky  and  Bertoni’s 
report  [2],  and  from  several  other  sources.  Where  possible,  we  chose  the  grid- 
point  values  so  as  to  preserve  the  large-scale  gradients  of  height. 

In  Figure  2,  the  distribution  of  the  surface -roughness  coefficient  used  in  the 
model  experiments  is  given.  These  values  were  subjective  estimates  made 
with  the  help  of  some  characteristic  values  of  this  quantity  for  various  ground- 
cover  characteristics  presented  by  Kung  and  Lettau  [23].  These  values  and 
their  distribution  are  of  course  open  to  considerable  criticism.  One  is  not  dis¬ 
pleased,  therefore,  by  the  result  of  the  subsequently  reported  experiment,  which 
indicates  that  this  coefficient’s  variation  may  not  be  overly  significant  for  large- 
scale  boundary-layer  processes. 
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VII 


EXPERIMENTAL  RESULTS 

25.  The  Synoptic  Case  and  the  Experiments  Conducted 

According  to  our  original  plans,  we  were  to  use  for  boundary  conditions 
at  the  top  of  the  boundary  layer  the  temperature,  relative  humidity,  and  geo¬ 
potential  fields  predicted  by  the  cloud-forecasting  technique  being  developed  by 

the  Computer  Techniques  Division  at  Headquarters,  Air  Weather  Service 
9 

(AWS)  .  We  were  notified  by  AWS  that  this  model  would  be  run  during  the  period, 
6  February— 11  February,  1964.  Consequently,  we  collected  synoptic  data  for 
this  period  and  determined  that  the  synoptic  case  studied  would  be  selected 
from  this  time  interval.  The  most  interesting  weather  phenomena  occurred  at 
the  beginning  of  this  period  in  the  region  indicated  in  Figure  1.  The  experiments 
discussed  in  this  section  were  all  twelve-hour  forecasts  for  the  period  1200Z, 

6  February— OOOOZ,  7  February,  1964. 

Synoptic  charts  for  the  initial  and  final  times  of  this  period  are  presented  in 
Figures  3  and  4.  Throughout  this  forecast  period,  low  cloudiness  covered  the 
greater  part  of  the  region  of  interest.  At  the  initial  time,  precipitation  was 
occurring  in  a  large,  crescent-shaped  area.  By  the  end  of  the  period,  the 
precipitation  was  confined  predominantly  within  the  northeastern  quadrant  of 
the  region,  with  a  southward-oriented  tongue  along  the  Appalachian  Mountains. 
This  development  was  associated  with  a  rapid  (25  knots)  eastward  displacement 
of  the  low-pressure  system  and  the  northeastward  movement  of  a  cyclogenetic 
area  originally  located  in  eastern  North  Carolina. 

The  basic  prediction  obtained  with  the  model  is  denoted  here  as  experiment 
number  two.  (An  error  in  the  specification  of  the  input  data  was  made  in  the 
first  experiment.)  Additional  experiments  were  made  with  the  modifications 


Refer  to  the  monthly  progress  report  for  June,  1964  for  the  reasons  for 
changing  this  plan.  (External  Monthly  Progress  Report,  Contract  No.  11910, 
The  Travelers  Research  Center,  Inc.). 
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Figure  3.  Synoptic  chart;  1200Z,  6  Feb.  1964.  - Sea  level  isobars 

700  mb  geopotential  contours  ti''/:'.?;:!  Broken  to  overcast  |  j  Precipitation 


indicated  in  Table  II.  All  the  experiments  utilized  thirteen  grid  points  in  the 
vertical.  They  were  located  at  50.  100,  150,  220,  300,  400,  500,  650,  850,  1150, 

1550  and  2000  m  above  the  terrain  at  the  grid  point.  A  total  of  one  hundred  grid- 
points  was  used  in  the  horizontal.  There  are  ten  columns,  numbered  from 
L  -  1  to  L  =  10,  from  west  to  east  and  ten  rows,  numbered  from  M  =  1  to 
M  =  10.  from  south  to  north. 

The  objective  of  these  experiments  was  to  estimate  the  general  predictive 
capability  of  the  model  and  to  determine  if  it  might  be  simplified  by  various 
modifications  without  affecting  its  predictive  accuracy.  In  what  follows,  we  dis¬ 
cuss  our  analvsis  of  these  experiments.  As  may  be  easily  visualized,  the  output 
from  such  a  series  of  experiments  is  voluminous.  The  approach  adopted  for  our 
analysis  of  results  is  basically  subjective.  It  must  be  clear  that  we  cannot  repro¬ 
duce  here  more  than  a  few  selected  figures.  We  hope  that  these  will  be  suggestive 
of  the  results  obtained  and  serve  to  illustrate  some  of  our  remarks. 

26.  The  Basic  Experiment  and  its  Verification 

Experiment  2  was  based  on  the  complete  model  outlined  in  the  preceding 
text.  The  only  significant  differences  between  the  model  outlined  there  and  that 
used  are  as  follows. 

First,  the  procedure  for  computing  the  various  eddy  fluxes  did  not  include 
the  free-convective  formulas.  Whenever  0^  was  computed  to  be  complex,  we  set 
the  radical  in  Eq.  IV-20  equal  to  unity.  This  should  have  the  effect  of  under¬ 
estimating  the  heat  flux  in  unstable  conditions. 

Secondly,  the  boundary  conditions  at  the  top  of  the  transition  layer  were 
taken  from  our  analysis  of  observations  at  the  initial  and  final  times  of  the  fore¬ 
cast  interval.  The  intermediate  values  were  specified  by  assuming  a  linear  vari¬ 
ation  in  time. 

The  prediction  made  in  experiment  two  is  indicated  in  Figure  5,  which  is  a 
constant-level  chart  showing  the  temperature  and  relative  humidity  forecast  at 
the  level  1500  m  above  mean  sea  level.  Figures  6  and  7  are  similar  charts  showing 
the  analysis  of  the  initial  and  verifying  data. 

Consider  first  Figures  5  and  6.  The  cloudiness  at  this  level  may  be  associated 
with  the  shaded  regions,  where  the  heavier  shaded  region  is  an  area  of  maximum 
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Figure  5.  Constant  level  analysis  (z 
OOOOZ,  7  Feb.  1964.  I  I  RH  <  77% 


=  1500  m);  Experiment  no.  2  forecast; 
RH  >  77%  RH  >  93% 


Figure  6.  Constant  level  analysis  (z  =  1500  m):  1200Z,  6  Feb.  1964. 

I  I  T-T^  >  3°C  CTl  T-T^  <  S-’C  1 . J  T-T^  ^  rC 
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Figure  7.  Constant  level  analysis  (z  =  1500  m);  OOOOZ,  7  Feb.  1964 
I  I  T-T^  >  3°C  mm  T-T  ,  <  3”C  I  I  T-T  ,  <  rC 


likelihood  of  overcast  cloud.  The  forecast  displacement  and  enlargement  of  this 
area  agrees  broadly  with  the  observed  change  in  cloudiness  and  precipitation  (see 
Figures  4  and  7).  The  temperature  forecast  in  the  central  portion  of  the  map 
is  also  commensurate  with  the  observed-temperature  variation.  Near  the 
boundaries  of  the  region,  the  forecast-temperature  field  is  in  considerable  error. 
This  is  discussed  in  more  detail  in  Subsection  31. 

Figures  8  and  9  depict  vertical  cross-section  analyses  along  the  sixth 
grid  row  (counting  from  one  at  the  southern-most  grid  row).  Again,  except  for 
the  region  near  the  left-  and  right-hand  boundaries,  the  temperature  and  humidity 
forecasts  agree  well  with  the  observed  distributions. 

In  Tables  III  and  IV,  the  predicted  vertical  velocities  are  given.  The  units 
are  cm  sec  and  positive  values  indicate  upward  moving  air.  The  terrain- 
induced  velocities  (w)  are  relatively  large  compared  with  the  frictionally-induced 
velocities  (w).  Note  that  it  is  the  large  value  of  w  at  L  =  7  and  M  =  6  which  is 
associated  with  the  column  of  maximum  shading  in  Figure  10.  Additionally,  we  may 
note  the  large  negative  value  of  w  at  grid  point,  L  =  8,  M  =  6.  The  importance 
of  this  terrain-induced  velocity  field  will  be  discussed  in  subsection  29, 

Table  V  presents  the  error  in  the  temperature  forecast  at  three  different 
levels  in  Experiment  2.  In  view  of  our  use  of  linear  interpolation  in  time 
between  observed  temperatures  as  the  upper-boundary  condition,  the  errors  in 
temperature  appear  to  be  quite  large.  Tables  VI,  VII,  and  VIII  present  the  grid- 
point  distribution  of  the  differences  between  observed  and  predicted  temperatures. 
The  largest  differences  are  noted  near  the  lateral  boundaries.  This  result  was 
expected  (see  Section  IV).  The  errors  and  their  distribution  may  be  accounted 
for  in  part  by  the  distribution  of  the  observations  used  to  specify  the  initial  and 
final  distributions  of  data.  In  addition,  the  linkage  between  the  upper-boundary 
conditions  (2000  m  above  ground)  and  the  interior  grid  points  is  predominantly 
associated  with  vertical  advection.  The  somewhat  greater  magnitude  of  the 
error  at  1550  m  indicates  that  this  process  is  not  overly  effective  in  forcing  the 
interior  temperatures  to  follow  those  applied  at  the  upper  boundary.  We  plan  to 
follow  up  this  implication  in  future  experiments  to  determine  the  possibility  of 
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Vertical  cross  section  (M  =  6);  Experiment  no,  2  forecast;  OOOOZ, 

IRH  ^  100%  ^  93%  fWPl  RH  ^  77%  RH  <  77% 


Figure  8 
7  Feb.  1964. 
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Vertical  cross  section  (M  =  6);  Experiment  no.  5  forecast;  OOOOZ, 
IrH  ^  100%  ^  93%  HHRH  >  77%  I  i  RH  <  77% 


TABLE  V 

DISTRIBUTION  OF  ERROR  (observed  less  forecast)  IN  TEMPERATURE  PREDICTED 
IN  EXPERIMENT  NO.  2  OVER  100  GRID  POINTS  AT  THREE  LEVELS  IN  THE 
VERTICAL  (measured  above  local  terrain),  AND  OBSERVED  TEMPERATURE 

CHANGES  (persistence) 


Class  Mean  (°C) 


Temperature 
error  (°C) 


Experiment  no.  2 


1550  m 


1150  m 


500  m 


1550  m 


Persistence 


1150  m 


500  m 


3.09 

3.02 

2.84 

4.40 

4.60 

3.90 
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dispensing  with  independent  forecasts  of  the  upper-boundary  temperature  and 
humidity. 

Tables  IX,  X,  and  XI  present  the  observed-temperature  changes  during 
the  forecast  interval  at  the  various  grid  points.  A  comparison  of  these  observed 
changes  with  the  errors  in  the  forecast  indicates  that  at  many  of  the  points  with 
large  observed-temperature  changes  the  predictions  are  remarkably  accurate. 

It  is  our  opinion  that  the  results  of  this  experimental  forecast  are  quite 
encouraging.  The  sjmoptic  case  used  did  not  lend  itself  to  a  detailed  evaluation 
of  the  potential  of  the  model  for  developing  low  cloudiness.  The  model  does 
seem  capable  of  translating  cloud  fields,  and  in  some  instances  (see  Subsection  29) 
of  dissipating  them. 

27.  Influence  of  the  Time  Step  Used 

In  our  analysis  of  the  computational  stability  of  the  model  in  Section  V, 
it  was  determined  that  linear  stability  was  assured  if  the  time  step,  At,  did  not 
exceed  the  time  required  for  an  air  parcel  to  traverse  one  grid  interval.  It  was 
indicated  that  the  time  step  could  be  as  long  as  forty  minutes,  if  the  grid  interval 
was  100  km  and  the  average  wind  40  m  sec  ^ .  A  question  remained  as  to  the 
accuracy  of  the  solution  which  would  result  if  the  time  step  was  made  relatively 
large. 

To  get  an  idea  of  this  influence,  we  compared  two  forecasts  which  differed 
only  in  the  time  step.  In  experiment  three.  At  was  set  to  thirty  minutes.  These 
predictions  were  compared  with  those  obtained  in  experiment  two  for  which  At 
was  fifteen  minutes.  The  largest  temperature  difference  between  the  two  fore¬ 
casts  was  0.6°C.  At  most  grid  points,  the  difference  was  ±  0.1°C.  The  humidity 
distributions  were  similarly  unaffected  by  this  change  in  time  step. 

The  conclusion  that  the  time  step  may  be  safely  increased  to  thirty  minutes 
must  be  examined  further  for  cases  in  which  turbulent  exchange  plays  a  greater 
role  than  in  the  S3moptic  case  studied  in  this  series  of  experiments.  It  does  seem 
safe,  however,  to  use  the  longer  time  step  in  those  instances  when  advective 
processes  predominate  in  the  boundary  layer. 
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28.  Influence  of  the  Roughness-coefficient  Variation 

The  roughness  coefficient,  z^,  was  assigned  a  constant  value  of  1.0  cm  in 
experiment  four.  The  principal  influence  of  this  modification  might  be  expected 
to  show  up  in  the  frictionally-induced  vertical-velocity  component,  w.  Table  XII 
gives  the  distribution  of  w.  Although  the  percent  changes  in  this  quantity  are 
relatively  large,  they  are  almost  negligible  in  comparison  with  the  difference 
between  w  and  w  (see  Tables  III  and  IV).  An  examination  of  the  temperature 
and  moisture  forecasts  showed  that  this  modification  had  negligible  impact  on 
the  predicted  values  of  temperature  and  relative  humidity. 

Tables  XIII  and  XIV  present  the  values  of  the  friction  velocity  (cm  sec 

-2  -1 

and  heat  flux  (millicals  cm  min  )  computed  at  OOZ,  7  February,  1964  in 
experiments  two  and  four.  Both  quantities  reflect  the  variation  of  between 
the  experiments.  The  fact  that  such  variations  were  not  effective  in  modifying 
the  predicted  fields  implies  that,  in  this  synoptic  case,  the  eddy- transport 
mechanisms  were  of  secondary  importance. 

29.  Influence  of  the  Terrain-induced  Vertical  Velocity 

In  experiment  five,  we  deleted  the  terrain-induced  vertical  velocity,  w, 
from  the  prediction  equations.  Figure  11  displays  the  predicted  temperature 
and  humidity  distribution  at  the  constant  level— 1500  m  above  mean  sea  level. 

It  is  to  be  expected,  of  course,  that  these  predictions  will  differ  from  those 
obtained  in  experiment  two  in  the  vicinity  of  the  major  terrain  features  (see  Fig¬ 
ure.  1).  That  this  is  true  is  shown  in  Table  XV,  which  gives  the  vertical  average 
of  the  absolute  value  of  difference  in  temperature  forecast  in  the  two  experiments. 

Figures  8,  9,  and  11  display  vertical  cross-section  analyses  of  the  temper 
ature  and  humidity  as  observed  and  forecast  in  these  two  experiments.  The 
rather  large  value  of  w  at  gridpoint  L  =  8,  M  =  6  in  Table  IV ;  was  associated 
with  the  dissipation  of  low  cloudiness  on  the  leeward  side  of  the  Appalachians. 
When  w  was  neglected  in  experiment  five,  high  humidities  were  predicted  in  this 
region.  It  seems  quite  clear  that  the  inclusion  of  this  process  of  terrain-induced 
vertical  motion  within  the  model  is  important. 
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Figure  11.  Constant  level  analysis  (z 
OOOOZ,  7  Feb.  1964,  |  |  RH  <  77% 


-  1500  m);  Experiment  no,  5  forecast; 
RH  >  77%  i:  1  RH  >  93% 
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30.  Influence  of  the  Variation  in  the  Mixing  Coefficient 

We  neglected  the  dependence  of  the  mixing  coefficient  upon  height  and 
stability  in  experiment  seven.  As  in  experiment  four,  this  modification  should 
influence  the  computed,  frictionally-induced  vertical  velocity,  w  (see  Table  XVI). 

The  remarks  made  in  subsection  28  may  be  reiterated  here.  Figure  12  displays 
the  predicted  temperature  and  humidity  fields  at  1500  m  above  mean  sea  level 
obtained  in  this  experiment.  There  are  only  minor  differences  between  this  fore¬ 
cast  and  that  obtained  in  experiment  two.  This  is  shown  more  quantitatively  in 
Table  XVII,  which  gives  the  vertically-averaged  absolute  value  of  the  difference 
in  temperatures  predicted  in  experiments  two  and  seven. 

The  implication  of  this  result  is  that  our  present  formulation  of  the  dependence 
of  the  mixing  coefficient  on  stability  and  height  is  not  reflecting  any  synoptically 
important  processes.  We  had  previously  examined  the  effect  of  this  dependence 
in  an  idealized,  two-dimensional  experiment,  with  the  same  result.  It  appears 
that  we  may  do  well  in  the  future  to  neglect  in  the  model  the  added  complication 
of  computing  this  dependence. 

31.  Influence  of  Lateral-boundary  Conditions 

Experiments  six,  eight,  and  nine  differed  from  experiment  two  in  only  the 
procedure  used  for  computing  the  horizontal  advection  on  the  lateral  boundaries. 

The  character  of  these  differences  is  outlined  in  Table  II. 

Tables  XVIII,  XIX,  and  XX,  give  the  vertical  average  of  the  absolute  value 
of  the  difference  in  temperature  predicted  in  experiment  two  and  in  these  other 
forecasts. 

In  both  Experiments  2  and  6,  we  computed  the  advection  parallel  to  the 
lateral  boundaries  while  neglecting,  in  certain  instances,  the  component  of  advection 
normal  to  the  boundary.  A  comparison  of  the  signs  of  the  errors  on  the  boundary 
indicates  that  Experiment  6  has  even  larger  errors  on  the.  eastern  boundary  than 
those  found  in  Experiment  2. 

In  Experiment  8,  we  neglected  both  components  of  advection  on  inflow 
boundaries.  It  appears  from  a  close  examination  of  the  boundary  errors  that  the. 
procedure  used  in  Experiment  8  would  have  improved  the  accuracy  of  the  fore- 
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Figure  12. 
OOOOZ,  7  Feb. 


Constant  level  analysis  (z  =  1500  m);  Experiment  no.  7  forecast; 
1964.  I  I  RH  <  77%  RH  >  77%  FTTl  RH  >  93% 
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cast  only  near  the  boundary  in  the  northeast.  On  the  other  boundaries,  the  results 
are  slightly  inferior  to  those  obtained  in  experiment  two.  Similar  comments 
apply  to  the  comparison  of  experiments  nine  and  two;  however,  the  complete  neglect 
of  advection  on  the  boundary  in  experiment  nine  led  to  noticeably  larger,  overall 
boundary  errors. 

In  spite  of  these  differences,  it  is  quite  clear  from  the  tabulated  differences 
that  the  boundary  influences  were  not  propagated  far  into  the  interior  of  the  region. 
None  of  the  four  methods  used  so  far  to  specify  lateral  boundary  advection  can  be 
highly  recommended.  It  will  be  necessary  to  analyze  alternative  schemes  as 
part  of  our  future  work. 
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APPKNDIX  I 


COQUDINATE  S\’STEM 

A  natural  cooi’dinate  system  with  which  to  represent  the  distribution  of 
atmospheric  i)rocesses  is  a  spherical  system  fixed  in  the  rotatint;  earth. 

Maurwit/,  |20|  derives  Newton's  second  law  of  motion  in  such  a  system.  He 
then  derives  tlie  simpler  form  of  these  equations  under  the  tanj;ent-planc  coordinate 
transformation. 

It  is  possible  to  simplify  the  equations  of  motion  in  spherical  coordinates 
by  the  nep;lect  of  certain  terms,  and  thereby  arrive  at  a  form  analogous  to  that 
achieved  by  the  tangent-plane  transformation.  We  refer  to  thi.s  system  as  the 
quasi-Cartesian  form  of  the  equations. 

In  order  to  make  clear  the  coordinate  system  used  in  this  paper,  we  dis¬ 
cuss  the  approximations  necessary  to  derive  the  quasi -Cartesian  form  of  the 
equations  from  their  original  form  in  spherical  coordinates.  We  then  introduce 
a  modified  spherical  coordinate  system  and  show  that  this  lead.s  to  a  modified 
quasi-Cartesian  form  of  the  equations. 

The  spherical  system  of  coordinates,  r'.  and  X,  has  its  origin  at  the 
center  of  the  earth.  The  plane  X  =  0  is  a  meridian  fixed  in  the  earth  containing 
the  axial  vector,  H,  which  represents  the  angailar  velocity  of  the  earth’s  rotation. 
The  plane,  0  =  0,  is  normal  to  Sl  and  coincides  with  the  equatorial  plane  of 
the  earth.  The  coordinate,  v',  denotes  the  radial  distance  of  the  point  from  the 
center  of  the  earth.  Let  r'  -  a  be  the  mean  radius  of  the  earth  and  therefore 
mean  sea  level.  Let  r'  =  E  (X,  0)  be  the  equation  representing  the  actual  disUince 
of  the  earth's  surface  teiu’ain  from  the  center  of  the  earth. 

Let  u.  V,  and  w'  be  the  components  of  linear  velocity  of  a  fluid  clement 
at  the  point  (r'.  0.  X).  Then 

,  dX 

u  -  r  cos  V  ~7: 

(It 

V  ..  (Al-1) 

dt 

dr' 
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The  particle  derivative  has  the  same  value  in  either  coordinate  system, 


dt  dt  dt  ’ 


and  by  definition, 


dG_^  9G  d0  8G  dr 

dt  9t  9A.  dt  ^  90  dt  ^  9r  dt  ’ 


From  Eqs.  AI-1  and  15  we  may  rewrite  Eq.  AI-20  as 


dG  _  9G  ^  u  9G  V  9G 
dt  9t  r'cos0  9X  r'  90 

It  is  convenient  to  define  w  as 

w  =  —  -  w'  — 


dt  9r  ’ 


(AI-19) 


(AI-20) 


(AI-21) 


(AI-22) 


Consider  next  the  transformation  of  the  expression  for  the  velocity  divergence 
given  by  Eq.  AI-5.  It  follows  from  Eq.  AI-22  that 


9w" 

9r' 


9w 


dE 


9r  9r'  |dt 

but  by  using  Eq.  AI-20,  this  becomes 
9w'  9w  9u  1  9E  9v  1 


9w  9u  _  _ 

9r'  9r  9r'  r'cos  0  9X 


9E  9y _ 

^  9r'  r' 


OE 

90 
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(AI-24) 


v_  ^ 
-2  90 


or 


9w^  _  ^  _9u_  1  9E  9v  9E  ^ 

9r'  9r  9r'  r'cos  0  9\  ^  r'  9r'  90  r'  dt  ’ 


(AI-25) 
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Using  Eq.  AI-17,  the  other  derivatives  in  Eq.  AI-5  can  be  evaluated. 


div  V  = 


3w 

ar 
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aE 
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aE 


r'  cos  4>  9^  9r'  r'cos  4> 

1  dv  av  1  aE 
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dv 
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dv 
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(AI-26) 


'  2w'  V  .  ,  1  dE 

+  — -  —  tan  (p  -  —  — . 

r  r  r  dt 

Noting  the  cancelation  effect,  Eq.  AI-26  reduces  to 

1 


-  aw 

div  V  =  —  +  ^ 

ar  r  cos  <p  ax 


au  1 

+ 

r 


dv  dE 

d(p^  r'  -  dt 


r  2w  V- 

+  I  — 7  -  ~  tan  (p 
[  r  r 


(AI-27) 


The  equations  of  motion  transform  easily,  of  course,  in  view  of  the  unique¬ 
ness  of  the  particle  derivative. 

To  recapitulate  then,  the  coordinate  system  used  in  the  text  is  essentially 
a  modified  spherical  coordinate  system.  In  the  various  equations,  we  have  neglected 
the  terms  involving  the  spherical  shape  of  the  earth  and  other  small  terms. 

We  employ  the  notation 


and 


dx  =  r  cos  <p  dX 
dy  =  r  d 
dz  =  dr 


u  =  r  cos  (p 


dt 


d0 

V  =  r  , 
dt 


w  = 


dr 

dt 


(AI-28) 


(AI-29) 


When  r  appears  undifferentiated  in  Eqs.  AI-28  and  29,  the  mean  value, 
r  =  a,  is  to  be  used.  We  point  out  here,  however,  the  assumption  introduced  in  the 
derivation  of  the  thermod3Tiamic  energy  equation  in  Subsection  7  of  the  text.  The 
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APPENDIX  II 


UPWIND  ADVECTION  APPROXIMATION 


1  ri 

We  define  X  (f)  ’  ,  ,  an  approximation  to  u  — ,  as  follows. 
■  'm,k  ax 


For  1  =  1; 


n  1,  n  li  n 

if  u  2;  0,  set  X  (f)  -  ,  =0, 
m,  k  m,  k 


ifu^’^,  <  0.  set  X  (f) 

m,  k  m,  k 


=  d  ^  u^’”  [f^'” 

m,k  m,k  |  m,k 


(AII-1) 


For  1  =  L  +  1; 


ifu^^^;”  s  0.  set  X(f)  '^^^’”  =  0, 
m,  k  m,  k 

L+l,n  L+l,n  .L+l,n  L,  n 

if  u  ,  >  0,  setX(f)  ,  =  d  u  ,  f  ,  -  f 

m,k  m,k  m,k  m,k  m,k 


(All- 2) 


For  1  <  1  <  L  +  1; 

if  u  ’  s:  0,  set  X  (f)  ’  =  d  u  ’  f  ’  -  f  ’  , 

m,k  m,k  m,k  m,k  ni,kj 

•  r  l.n  ^  1. Ffl+l.n  -l,n 

if  u  ’  <  0,  set  X  (f)  ’  =  d  u  ’  f  -  f  , 

m,k  m,k  m,k  m,k  m,  k 


(AII-3) 


We  define  Y(f)  ,  an  approximation  to  v  as  follows, 
m,  k  9y 


For  m  =  1; 


i£v;-_"  ^  0.=etY(f);;;  =  0. 


if  V  ”  <  0,  set 
1,  k 


Y(f)^’”  =  d  ^  v^’”  T^’”  -  f^’” 

^  ^l,k  l,k  !  2,k  l,k 


(AII-4) 


For  m  =  M  +  1; 

ifvj’”  <  0,  setY(f)i’”  =0, 

M+l,k  '  'M+l,k 

^  '^M+l.k  ^M+l,k  "  ^M,k 


(AII-5) 
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For  1  m  M  +  1; 


In  l,n  ,-l  l,n 

if  V  ’  >  0.  set  Y(f)  ’  ,  =  d  V 

m,k  ^  ^m,k  m,k 

<  0,  set  =  d  ^  ! 

m,k  ’  ^^m,k  m,k 


l,n 


rl>n 

m-1,  k 
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^l,n  _  jl,n 
m+l,k  m,k 
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APPENDIX  III 

THE  DIFFERENCE  APPROXIMATION  FOR  THE  HEAT-DIFFUSION  EQUATION 


Set  t  =  t,  x  =  x,  ,  y  =  y. 


For  k  =  1; 


'  u,l 


X(T)^’  J'  -  Y  (T)^’''  -  f  - 


(AIII-l) 


X,  i^+1  _  )^+l 

_  _fc2,  ^  _E_  _ 

^  .  11,2  [  (Az)^  CpJ  V  v 


For  k  =  2,  ...,  K  -  1; 


T^,  <^+1  _ 

u,k  a,  k 
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APPENDIX  IV 


THE  DIFFERENCE  APPROXIMATION  FOR  THE  VAPOR-DIFFUSION  EQUATION 


Set  t  =  t,  y  =  y  ,  x  =  x,. 

u’  ^  X 


For  k  =  1; 


X,  j^+1  X,  u 


{-  X,  P+l  _  X,  jx+l- 
M,2  (Az)^  J 


For  k  =  2,  ...,  K  -  1; 
X,  i^+l  X,  V 


a  W  r 

+  (Az)j^  ["^M.k+l 


(AIV-1) 


X,  u+1  X,  1^+1 
H,k+1  ^iJL,k-l 


t  -  X,  u+1  _  X,  u+1 

^  _ 2  u,k+l  ^M.k 

(Az)k^l  +  (Az)k  \  M,k+1  (Az)k^j 


r  X,u  X,  u+1 
.X,^  PM,k  ^M,k-1 


(AIV-2) 
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APPENDIX  VI 


DEFINITION  OF  COEFFICIENTS  IN  EQUATION  V-3 


For  k  =  1; 

X, 

a  -  =  0 

M.  1 


=  -  {(Atf^  +  2X^*2  /  [(Az)^]} 
fu,  1  M>  z  z 

=  2K^’"  /  [(Az)^] 

M.  1  M.  2  2  . 


(AVI-1) 


X.j.  ^  ^  X(T)^’''  +  Y(T)^’''  +  ^  -  2gK^’^  /  [C  (Az)] 
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